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Chapter 1 



Introduction 



1.1 Introduction to the equations 

One of the classic problems of physics is that of the motion of a particle. 

In finite dimensions, i.e., in M", we have a particle described by its 
position X and velocity x, and the particle is under the influence of a poten- 
tial V{x) (corresponding perhaps to gravity or an electrical potential). The 
Lagrangian is the kinetic (energy of motion) minus the potential and is given 

by 

L — -m\x\'^ — V{x). 
The Euler- Lagrange equations give the equations of motion which are 

mx — — W{x). 

In field theory, the position of the particle is not given by a vector 
a; e M" but rather by a function which is a point in an infinite dimensional 
vector space. Call this position f{x) and instead of the inner product on 
R", the inner product is based on a function space, such as the £^ inner 
product. Using the £^ norm as an example, if there is a potential V{f ), and 
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the Lagrangian action is given by: 

\f{x)\'-V{f{x))dx, 



the calculus of variations gives the equations of motion 

Sf 

We will study two situations similar to this in this dissertation. 

First, we address the situation of the Yang Mills Lagrangian in 4 di- 
mensions. The Yang Mills equation is a generalization of Maxwell's equations 
in a vacuum. We wish our particles to have certain internal and external 
symmetries, which give rise to the various geometrical objects in the prob- 
lem. The state of our particles are given by gauge potentials or connections 
A on R^, where we identify as H, the quarternions. The gauge potentials 
have values in the Lie algebra of SU(2) which can be viewed as pure imagi- 
nary quarternions, Im(EI). The curvature Fij = diAj — djAi + [Ai,Aj] where 
[Ai, Aj] = AiAj — AjAi is the bracket in the Lie Algebra, gives rise to the 
potential V{A) = {F, F) which is a nonlinear function of A. The action is: 



L = IJ {F,j,F,,)dx. (LI) 



The local minima of [|L1|| are the instantons on 4 dimensional space. These 



correspond to solutions of Maxwell's equations in the vacuum, if we consider 
the Yang Mills Lagrangian to be a generalization of Maxwell's equations. We 
now consider the wave equation generated by this potential with Lagrangian: 

L = \j {A,, A,) - ^{F,j,F,,)dx (L2) 



Here we use the summation convention; however in taking inner products of 
two forms, we only sum over i < j, and not all i,j, which has here been 
expressed by dividing by 2 in computing {Fij,Fij). Now use the calculus of 
variations by taking A ^ A + 6 A, to obtain 

L{A + SA) = 

J ^{Ai, Ai) + {Ai, sAi) + ^{sAi, sAi) - 

^{diAj + diSAj - djAi - djSAi + [Ai + 5Ai, Aj + 5Aj], 
diAj + diSAj - djAi - djSAi + [Ai + 6Ai, Aj + SAj])dx 
Using the definition of F, one obtains 

L{A + 5A) ^ J A,) + {Ai, SA,) + ^{6Ai, SAi) - ^{F,„ Fij) 

+^{F,„6F,,) + ^{6F,„6F,,)dx. 
One now takes the linear part of the variation, 

J {Ai,SAi) + ^{Fij,SFij)dx, 
and analyzes it. First work on 

J {Fi„ SF,,) = J {Fi„ diSAj) - {F,„ d,SA,) + (F,,-, [SA, A,]) + {F,„ [A„ SA,]) 

= J{F,„d^SAj)-{F,j,djSA,)- 

Tr{Fij{SAiAj - AjSAi + AiSAj - SAjAi)) 
= J -2{djFij, SAi) - TtiFijiSAiAj - AjSAi)) 

by integration by parts and the cyclic property of traces: 



4 

= J -2{djFij, 6Ai) - TriFi^A, - FijAj))6Ai 

= J -2{V,F,„6A,). 

Here Vj represents the covariant derivative in the j direction. Now integrate 
the first term by parts to obtain: 

J {Ai, SAi)dx = - J {Ai, 5Ai)dx. 

Put these two things together to get 

Ai = -V,-F,,- (1.3) 

This is our evolution equation for the 4+1 dimensional model. 

The second situation to be addressed is that of the CP^ model in 2+1 
dimensions. 

In the finite dimensional analog of this case, we consider motion on a 
manifold under the influence of a potential. The Lagrangian in this case is 

L = ]^gi3{x)xiXj - V{x), 

where Qij is the metric tensor that allows us to measure length on the curved 
space of the manifold. The equations of motion are of the form 

where is a covariant derivative with respect to time. 
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In the CP^ model we have motion on an infinite dimensional manifold 
as we consider maps from M^"*"^ ^ 5*^, to the two sphere, of a particular degree. 
If (/) : ^ ^2 ^jjg Lagrangian is 

Identifying S'^ as C + {oo}, one can consider u : M^"*"^ C + {oo}. The 
Lagrangian is 





\u 


2 






u\ 


2 


(1+1 


u\ 


(1 + 


\u 


2^2 



(1.4) 



One can see this is one again the kinetic minus the potential energy. 

The calculus of variations on this Lagrangian in conjunction with inte- 
gration by parts yields the following equation of motion for the CP^ model: 

(1 + \u\''){d^u - dlu - dlu) = 2u{\dtu\^ - \d^uf - \dyu\^) (L5) 
1.2 Introduction to the adiabatic limit 

The origin of the term "adiabatic" seems to be with reversible processes: pro- 
cesses that do not change the entropy of a system and consequently that can be 
undone. One way of enacting such a reversible process it to make the changes 
infinitely slowly. 

Another way of looking at making changes infinitely slowly is to break 
space-time into two pieces, space and time. If one then rescales the spatial 
piece to be extremely small, or the time piece to be extremely large, it then 
takes correspondingly more time to move from point to point. In the limit as 
we take the spatial piece to be infinitely small or the time piece to be infinitely 
large, we are moving infinitely slowly. 



Naturally, this sort of rescaling can occur in any problem that can be 
broken into two pieces. Any space that is a Cartesian product, say an infinite 
cylinder, 5"^ x R or a torus x S^. Or even any space that locally looks like 
a product, such as our old friend the Mcebius band, or, if one is even more 
adventurous and brave, a twisted bundle. The main requirement is that we 
have a metric on the space that splits into two pieces, 

g(x,y){x,x) + h(x,y){y,y). 

Then one can introduce a parameter A without changing the topology of the 
space: 

g{x,y){x,x) + Xh{x,y){y,y). 

An adiabatic limit is simply taking the parameter A ^ or A ^ oo, which 
makes the space Y look either infinitely large or infinitely small. While chang- 
ing from one finite and nonzero A to another causes no changes in the topo- 
logical properties of the space, the limit has no such immunity. Conclusions 
based on an adiabatic limit must be examined with care, for what happens in 
the limit might not be close to what happens near the limit. 

Adiabatic limits take on special meaning in the case of space-time. 
Taking an adiabatic limit results in moving either infinitely slowly or infinitely 
quickly, usually infinitely slowly. When we have a partial differential equation 
with known moduli spaces of static solutions, the adiabatic limit tells us that 
as velocity tends towards zero, we should get motion along the geodesies of 
these moduli spaces. We approximate solutions with small velocity with these 
geodesies, and hence this is called the Geodesic Approximation. Investigating 



7 

this phenomenon in the two cases mentioned in the introduction is the concern 
of this dissertation. 



Chapter 2 



The 4+1 dimensional model 



The first things to identify in this problem are the static solutions to equation 



||1.3|| . These are simply the 4 dimensional instantons investigated in [|l|. In 
the form of all such instantons in the degree one sector is shown to be: 

, , , 1 ( (x — a)dx — dxix — a) 1 , 

A{x) = - { ^ — \- — } x = xi + x2i + xsj + Xik e m. 

2 [ X'^ + \x — J 

The curvature F of this potential is computed by F = dA + [A, A] and is: 

dx Adx 
^ (A2+ |x-a|2)2- 

One notices the denominators are radially symmetric about x = a. 

An instanton of unit size centered at the origin would be 



A{x) 



|2 



1 + \x\ 

One motion one can study from these static instantons would be to consider 



connections the form 



1 I ^ijdiOc di^jCijC 

2 l/(r,t)+r 



A{r, t) = - I \ r = \ xl + xl + xl + x\ 



and derive an equation of motion for f{r,t) from ||1.3|| . This is actually quite 
computationally extensive. 
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To get at this, start with connections of the form 



A{r, t) = -^9^^' {^^^ ~ dxx} . 

Now we have a formula for F given hj F = dA + [A, Al\. But F cannot easily 
be expressed as before. The components can be computed using Maple, with 
quarternions defined as 



i 

1 



3 



-1 

1 



k 



i 
-i 



the usual derivative, and the bracket is that of two matrices: [M, A^] = MN ■ 
NM. First, computing F one obtains components such as 

'xlg'{r, t) + 2r^(r, t) + xlg\r, t) 



'12 



2xig{r,tY-2xig{r,ty 



+ 



r 

g'{r,t) 



+ 2g{r,tfj (X2X3 - X1X4) 
+ 2g{r,tf ) (xia;3 + a;2a;4) 



+ 



k 



Once F has been obtained one may apply ||1.3|| to obtain an equation for g{r, t). 
After much labor, one obtains: 



5a' 

9 = I2g^ + — + 9"- 8g^r^. 

r 



:2.i] 



if 



It is 



It is much less difficult to now compute a differential equation for /(r, t) 



9{r,t) 



/(r,t) + r2 



The static solutions for f{r,t) are simply horizontal lines, f{r,t) = c. 
The adiabatic limit expects that motion under small velocities should progress 
from line to line, i.e., f{r,t) = c{t). f{r,t) = is a singularity of the 
system, where the instantons are not well defined. We can use this numerical 
approximation to the adiabatic limit to observe progression from /(r, t) = cq > 
towards this singularity. Our initial assumption is that f{-,t) = c{t). 



2.1 Numerics for the 4+1 dimensional model 



A finite difference method is used to compute the evolution of ||2.2|] numerically. 



Unless otherwise noted, centered differences are used consistently, so that 

g{x + 5) - g{x-5) 



9\x) 
9"{x) 



25 

g{x + 5) +g{x -5)- 2g{x) 



In order to avoid serious instabilities in 2.21, the terms 



f" + — (2.3) 
r 



must be modeled in a special way. Allow 



r+ — = 

r 



where 



In Appendix A, we will see that this operator has negative real spectrum, 
hence it is stable. However the naive central differencing scheme on [ ^.3 



always results in uncontrolled growth near the origin. General wisdom holds 
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that when one has difficulties with the numerics in one part of a problem one 
should find a differencing scheme for that specific part in the natural to that 
specific part. Applying this allowed for the removal of the problem near the 
origin. Instead of using centered differences on /" and on /', we difference the 
operator: 



The "natural differencing scheme" is 



Lf ^ r" 



^ I nV /(r + ^)-/(r) 



5 

r 

2 



/(r) - /(r - 6) 
6 



With the differencing explained, to derive f{r,t + At), one always has 
a guess for /(r, t + At) given by either the initial velocity, e.g. /(r, t + At) = 
/(r, t) + voAt, or by /(r, t + At) = 2/(r, t) - f{r, t-At). Use this to compute 
f{r,t) on the right hand side of ||2.2|| . Then solve for f{r,t + At) in the 
difference for f{r,t), and iterate this procedure to get a more precise answer. 
So, one iterates 

5/'(r,t) 



f{r,t + At) = 2/(r,t)-/(r,t- At) + (At)2 

2/(r,t)2 2r(r,t)2 



r(r,t) 

8f'{r,t)r 



f{r,t)+r'^ f{r,t)+r'^ f(r,t) + r'^ 

where all derivatives on the right hand side are represented by the appropriate 
differences. 

There remains the question of the boundary conditions. The function 
/ is only modeled out to a value r = R ^ 0. Initial data for /(r, 0) was 
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originally a horizontal line. The corresponding boundary conditions are that 
/(i?, t) = f{R- Ar, t), and that 

/(0,t) = ^/(Ar,t)-i/(2Ar,t) 

i.e., that / is an even function. 

Subsequent investigation of the model indicated that the appropriate 
form for f{r,t) was a parabola instead of a line, and the f{R,t) boundary 
condition was changed to reflect this. For the runs with parabolic initial data, 
we set the boundary condition at i? to be: 

R 



f{R,t)=f'{R-Ar,t) 



R- Ar' 



2.2 Predictions 



Equation ||1.2|| gives us the Lagrangian for the general version of this problem. 



1 [ ^jcdi'jc di^Xjtic 



We are using 



and the adiabatic limit says we will move on the moduli space of these solu- 
tions. This will give us an effective Lagrangian. The portion of the integral 
given by 



1 

represents the potential energy and integrates to a topological constant, hence 
it may be ignored. We need to calculate 



{Ai, Ai)dx. 
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First calculate 

3^2 p 

So the effective Lagrangian is 



/ 



Letting y = x/ \ff we can rewrite this integral as 

3/' S \V? 



f - 

(1 



The integral with respect to y converges, hence we have 

This is purely kinetic energy. Since the potential energy is constant, so is the 
kinetic energy. We have 

7 = '- 

Integrating this we get 

If / = occurs at time T, we find 

T = - — 

Cl 

hence we rewrite this as 

f^a{t-Tf. 



This is how we predict that f{0,t) will evolve. 
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2.3 Results 

The computer model was run under the condition that /(r, 0) = /o with various 
small velocities. The initial velocity is /(r, 0) = fo, other input parameters 
are R = Vmax, Ar and At. 

2.3.1 Evolution of /(O, t) 

The first question to ask is how does the evolution of the origin occur. We 



note that equation ||2.2|| becomes the regular linear wave equation 

/ = /" 

as r ^ CX3, and so the interesting nonlinear behavior is at the origin. Conse- 
quently we track the evolution of /(0,t). 

The evolution of /(0,t) is a parabola of the form: 

f{0,t)=a{t-Tf, 

where 

4/„ 



and 

T = — 

1^0 r 

A typical evolution of f{0,t) is given in Figure In this figure, the 
equation 0.000025(t — 200)^ neatly overlays the graph of f{0,t). This picture 
represents the evolution where /o = 1.0 and vq = —0.01. Hence c — ^MlL — 



4(1.0) 

0.01 



0.000025 and T = = 200, as predicted 



The time to "blow up" is the parameter T in this equation. Recall that 
fo = f{r, 0) is the initial height and /o(r, 0) = vq is the initial velocity. Using 
a least squares parabolic fit to the origin data obtained after f{0,t) < 0.5/o, 
one obtains the parameters a and T for a given origin curve. Table |2.1| shows 
the behavior. 
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Evolution of the Origin 




-0.05 I 1 1 1 1 1 1 1 1 

40 60 80 100 120 140 160 180 200 



Figure 2.1: 4+1 dimensional model, evolution f{0,t). 



Table 2.1: 4+1 Dimensional Model, parabolic fit to /(O, t) vs. Initial conditions 
/o and vo 



fo 


^^0 


a 


T 


1.0 


-0.010 


0.00002501 


200.1 


2.0 


-0.010 


0.00001257 


399.4 


0.5 


-0.010 


0.00005166 


99.0 


4.0 


-0.010 


0.00000626 


799.7 


4.0 


-0.020 


0.00002503 


400.1 


4.0 


-0.005 


0.00000157 


1599.3 
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2.3.2 Characterization of time slices /(r, T): evolution of a horizon- 
tal line 

The most striking immediate result is that the initial line, /(r, 0) = /o, evolved 
an elliptical bump at the origin that grew as time passed. Figure ^]2| shows 
this behavior. 

The elliptical bumps can be modeled as 

{y - kY 



a 



The question naturally arises as to how the parameters a, b and k evolve. This 
is straightforward: 

a = t 

4/o 
k = fo + vot 

Table shows these values as calculated using least squares fitting (to either 
a line or parabola). Since the elliptical fit only works before the right end of the 
ellipse hits the boundary at r = i? the fit sometimes needed to be restricted to 
the portion of the data before this occurred. While the ellipse is small, there 
is a great deal of noise in finding the elliptical parameters, and to get a good 
data fit, this noise must often be removed. Here and ba is the slope and 
intercept of the line a{t), and likewise and bk are the slope and intercept 
of the line k{t). The parameter c is that in b{t) = ct^ 

A typical evolution for a with /o = 4.0 and vq = —0.01 is in Figure p.3| . 
Note the noise when t < 100. 



Figure 2.4 is a typical evolution for b with /o = 4.0 and vq = —0.02. 
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Figure p]5| is a typical evolution for k with /q = 2.0 and Vq = —0.01. 
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0.8 
0.6 
0.4 
0.2 
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-0.4 
-0.6 h 
-0.8 



Time Profiles of Data 



T = 20s 
T = 40s 
T = 60s 
T-80S 
T = 1 00 s 
T = 12Qs 
T = 140 s 
T = 1 60 s 
T = 180 s 



20 40 60 80 



100 120 140 160 180 200 
r 



Figure 2.2: 4 + 1 Dimensional model, Time Slices /(r, T) evolve an elliptical 
bump at the origin 



Tabic 2.2: 4+1 Dimensional Model, Elliptical Parameters vs. Initial Condi- 
tions /o and vq 



k 


^^0 


rria 


ba 


c 


rrik 


bk 


1.0 


-0.010 


0.999 


-0.11 


0.00002510 


0.0100 


-1.00 


2.0 


-0.010 


0.999 


0.01 


0.00001270 


0.0108 


-2.00 


0.5 


-0.010 


1.007 


-0.420 


0.00005653 


0.0100 


-0.51 


4.0 


-0.010 


1.001 


0.764 


0.00000626 


0.0100 


-4.001 


4.0 


-0.020 


0.998 


0.177 


0.00002527 


0.0200 


-4.003 


4.0 


-0.005 


1.001 


-0.834 


0.00000157 


0.0050 


-4.000 




Figure 2.3: 4+1 dimensional model: elliptical parameter a as a function of 
time when Jq = 4.0 and Vq = —0.01. 




Figure 2.4: 4+1 dimensional model: elliptical parameter 6 as a function of 
time when Jq = 4.0 and Vq = —0.02. 
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Elliptical parameter k as a function of time 



k(t) 




50 100 150 200 250 300 350 400 

time in seconds 



Figure 2.5: 4+1 dimensional model: elliptical parameter /c as a function of 
time when /q = 2.0 and Vq = —0.01. 
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2.3.3 Characterization of time slices /(r, T): evolution of a parabola 

The elliptical bump that formed in the evolution of a horizontal line and the 
various configurations that ensued after it bounced off the r = and r = R 
boundary suggested that perhaps the curve was trying to obtain the shape 
of a parabola. After all, near r = 0, ellipses are excellent approximations for 
parabolas of the form 

f{r,t)=pr^ + h. (2.5) 
To get the parabola, calculate from the general form of our ellipse in 

21 

dy x^lP' 
dx {y — k)a? 

so 

d'^y —6^ xh"^ dy 

dx {y — k)a'^ {y — kya^ dx 

At x = 0, y — k = b and this gives 

d'^y —b 
dx 0? 

Recall from the previous section that b = ct^ and a = t, so this gives 

dx'^ 



The identification of c gives 



<Py_ ^ 

dx^ 4/o 



So 



When a run is started with this initial data, fo = Vq = —0.01, /o = 

2 

/(O, 0) = 1.0 and p = = —0.0000125, the time shces of the data have this 
same profile. This is shown in figure pTB. 



The curvature of the parabola at the origin, as measured by the pa- 
rameter p from equation ||2.5|| changes by less than 1 part in 100 during the 



course of this evolution, a graph of p over time can be seen in figure pTf . 

The other parameter in equation [ ^.5| for the parabola, h, should be 
given by the height of the origin, but this was calculated in the previous section 
to be a(t — T)^, substituting the expressions for c and T we obtain: 



4/o V 1^0 1 

This is indeed the correct form, as shown in figure |2.8| . The initial conditions 
were Vq = -0.01 and fo = /(O, 0) = 1.0, hence h{t) = 0.000025(t - 200)^. The 
plot of the function overlays the data. 

Now, using both expressions for p and h, one can get the general form 
of a parabolic f{r,t), which is 

^0 2 , ^0 f. 2/o 



8/o 4/o V mJ 
Substitute this into the partial differential equation [p.2|| , get a common de- 
nominator and simplify to obtain: 
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Clearly these two sides are not the same, and the error between them is 



Since our concern is the geodesic approximation, I'o/l/o) is always chosen to 
be less than 1/200. This term is always much smaller than the 2r^ term, 
and since we only have accurate numerics away from the neighborhood of the 
singularity, the term with t — 2/o/|vo| is on the order of magnitude /o/lfol, 
which is large compared to the correction. 
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Time profiles with paraboiic initial data 



t = 20s 

t = 40s 



t = 80s 

= 100 s - - - 

= 120s - - - 

= 140s 

= 160s 

= iso-s-^— - 



10 



20 



30 



40 



50 



Figure 2.6: 4 + 1 dimensional model: Time slices of the evolution of a parabola 
are parabolas. 
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Parabolic parameter p as a function of time 




40 50 60 
Time in seconds 



100 



Figure 2.7: 4 + 1 dimensional model: evolution of parabolic parameter p with 
time. 




Figure 2.8: 4+1 dimensional model: Parabolic parameter h as a function of 
time evolves as f{0,t). 



Chapter 3 
The CP^ model, charge 1 sector 



The first thing to identify in this problem are the static solutions determined 
by equation ||1.5|| . These are outlined in among others. The entire space of 
static solutions can be broken into finite dimensional manifolds M„ consisting 
of the harmonic maps of degree n. If is a positive integer, then consists 
of the set of all rational functions oi z = x + iy oi degree n. For this chapter, 
we restrict our attention to Mi, the charge one sector, on which all static 
solutions have the form 

u = a + I3{z + ^)-\ (3.1) 
In order to simplify, consider only solutions of the form 

or look for a real radially symmetric function f{r,t) so that this evolution 
occurs as: 

f{r,t) 

z 

It is straightforward to calculate the evolution equation for f{r,t). It 

is: 
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The static solutions for f{r,t) are the horizontal lines f{r,t) = c. In 
the adiabatic limit motion under small velocities should progress from line to 
line, i.e f{r,t) = c(t). f{r,t) = is a singularity of this system, where the 
instantons are not well defined. We use this to form a numerical approximation 
to the adiabatic limit to observe progression from /(r, 0) = cq > towards 
this singularity. 



3.1 Numerics for the CP^ charge 1 sector model 

A finite difference method is used to compute the evolution of ||3.2|] numerically. 
As with the 4 + 1 dimensional model, centered differences are used consistently 
except for 

f" + — . (3.3) 

r 



In order to avoid serious instabilities in ||3.2|| this is modeled in a special way. 
Let 



Lf = r-'dydrf = r 



r 



In Appendix B it is shown that this operator has negative real spectrum, hence 



it is stable. The naive central differencing scheme on ||3.3|| results in unbounded 
growth at the origin, but the natural differencing scheme for this operator does 
not. It is 



31 



With the differencing explained, we want to derive f{r,t + At). We 
always have an initial guess at /(r, t + At). In the first time step it is /(r, t + 
At) = f{r,t) + foAt with vq the initial velocity given in the problem. On 
subsequent time steps /(r, t + At) = 2/(r, t) — /(r, t — At). This can be 
used to compute /(r, t) on the right hand side of ||3.2|| . Then solve for a new 
and improved /(r, t + At) in the differencing for the second derivative f{r,t) 
and iterate this procedure several times to get increasingly accurate values of 
/(r,t). 

There remains the question of boundary conditions. At the origin /(r, t) 
is presumed to be an even function, and this gives 

/(0,t) = ^/(Ar,t)-^/(2Ar,t). 

At the r = R boundary we presume that the function is horizontal so f{R, t) = 
f{R-Ar,t). 

3.2 Predictions 

Equation ||1.4|| gives us the Lagrangian for the general version of this problem. 



We are using 

u = — 

z 

for our evolution, and via the geodesic approximation, we restrict the La- 
grangian integral to this space, to give an effective Lagrangian. The integral 
of the spatial derivatives of u gives a constant, the Bogomol'nyi bound, and 
hence can be ignored. If one integrates the kinetic term over the entire plane, 
one sees it diverges logarithmically, so if /3 is a function of time, the soliton 
has infinite energy. 
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Nonetheless, this is what we wish to investigate. We cannot address 
the entire plane in our numerical procedure either, hence we presume that 
the evolution takes place in a ball around the origin of size R. II (3 — f{r, t) 
shrinks to in time T, we need R > T. Under these assumptions, up to a 
multiplicative constant, the effective Lagrangian becomes 



/ rdr- 
Jo I 



,2f2 



2\2 



which integrates to 



In 1 



i?2 



Since the potential energy is constant, so is the purely kinetic Lagrangian, and 



2 



In 1 + 



R^ 



R^ 



2' 



with c (and hence c^/2) a constant. Solving for / we obtain 



R^ 



(3.4) 



P + R\ 

Since we are starting at some value /o and evolving toward the singularity at 
/ = this gives: 



R^ 



cdt. 



(3.5) 



The integral on the right gives ct. The integral on the left can be evaluated 
numerically for given values of R, fo and /(O, t). A plot can then be generated 
for ct vs. f{0,t). What we really are concerned with is f{0,t) vs. t, but once 
the value of c is determined this can be easily obtained. One such plot with 
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/o = 1.0, R = 100 of /(O, t) vs ct is given in Figure p.l| . This curve is almost, 
but not quite, linear, as seen by comparison with the best fit hne to this data 
which is also plotted in Figure The best fit line is obtained by a least 
squares method. 
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3.3 Results 

The computer model was run under the condition that /(r, 0) = /o with various 
small velocities. The initial velocity is /(r, 0) = fo, other input parameters 
are R = Vmax, Ar and At. 

3.3.1 Evolution of /(O, t) 

The primary concern with the evolution of the horizontal line is the way in 
which the singularity at /(0,t) = is approached, because once again as 
r ^ oo equation ||3.2|| reduces to the linear wave equation 

/ = /", 

and so we expect the interesting behavior to occur near r = 0. The model 

is run with initial conditions that /(r, 0) = /o, f{r,0) = Vq. Other input 
parameters are R = r^^x, Ar and At. 

The evolution of the initial horizontal line seems to remain largely flat 
and horizontal, although there is some slope downward as time increases. This 
is shown in figure [3^ . 

We track /(0,t) as it heads toward this singularity, and find that its 
trajectory is not quite linear, as seen in figure ^]3[ This is suggestive of the 
result obtained in the predictions for this model. 

A simple characterization this evolution can be obtained via the slope 
m of the best fit line (via least squares) and the time T at which /(O, T) = 0. 
Step sizes of Ar = 0.01 and At = 0.001 were used with R = 100.0. Table ^ 
displays the initial data /o = /(r, 0) and Vq = f{r, 0) along with the consequent 
T and m values. 
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We have 



and 



T 



m 



l-2/o 

1^0 1 



A better method would be to try to use the resuh from equation p3 



in chapter 3.2. This requires a determination of the parameters R and c. We 
aheady have /o and /(O, t). To determine R and c, observe from equation 
that: 



1 



In 1 + 



R^ 



R' 



/2 + i?2 



(3.6) 



Since R is large and / is small 



In 1 + 



and 



R^ 



In 



/2 + i?2 



1. 



Consequently we can rewrite equation ||3.6|| as 



[ln(i?2) _ ln(/2) _ 1] 



The plot of ln(/) = ln(/(0,t)) vs 1/p = l/P{Q,t) should be linear with the 
slope m = and the intercept h = (21n(i?) — l)/c^. Such a plot is easily 
obtained from the model, and given slope and intercept, the parameters c and 
R are easily obtained. 

Figure is a plot of ln(/(0,t)) vs. l/p{0,t), with initial conditions 
Ar = 0.01, At = 0.001, /o = 1.0 and vq = -0.01. It is easily seen that 
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although the plot of ln(/(0,t)) vs. l/p{0,t) is nearly straight, it is not quite 
a straight line. This indicates that the values of R and c are changing with 
time. 

The best fit line y = mx + b has slope m = —2810 and b = 10200. We 

have 

V m 

and 

R = exp ( + i ) . 

This gives values c = 0.0267 and R = 62.1. Using these values of c and R in the 
calculation of equation | p.5| |, we obtain the plot of f{0,t) vs t given in Figure 
|3.5| . This is overlayed with the model data for /(0,t) vs. t for comparison. 
These two are virtually identical. 

The next question is how do the parameters c and R vary with the initial 
conditions. The first piece of bad news is that they are strongly dependent on 
the size of Ar when Ar > 0.01, however, they appear to converge as Ar —>■ 0. 



Table 3.2 contains data of c, R vs. Ar and At under initial conditions /o = 1.0 



and fo = —0.01. 



Table |3]^ contains the data for c and R vs. change in the initial velocity 
Vq, under the initial conditions /o = 1.0, Ar = 0.01 and At = 0.001. The data 
for R varying with vq fits well to the parabola 

-0.00237878x2 + 0.73551687x + 3.23010905. 



This fit is shown in figure |3.6| . We know this must break down for small initial 
velocities, because it would send the value of R to — oo. The values of R do 
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seem to be levelling out. This suggests a direction for further research. 



Table |3]J containing the data for c and R vs. change in /o, the initial 
height. The parameter R varies close to linearly with /o, while the parameter 
c remains nearly constant, changing by less than 7% over the course of the 
runs. 
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Time Profiles of Data 



T = 1 2 sec 

T = 24 sec 



T = 36 sec 
T = 48 sec 
T = 60 sec 
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T = 84 sec 
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T= 120 sec 



0.8 



0.6 - 
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Figure 3.2: CP^ model, charge 1 sector: Time slices f{r,T) 



Table 3.1: CP^ model, charge 1 sector: Parameters for best fit line to f{0,t) 
vs. initial data /o and Vq 



fo 


^^0 


T 


m 


1.0 


-0.010 


113 


-0.00815 


2.0 


-0.010 


228 


-0.00811 


3.0 


-0.010 


346 


-0.00793 


4.0 


-0.010 


466 


-0.00836 


5.0 


-0.010 


588 


-0.00828 


1.0 


-0.020 


58 


-0.01586 


1.0 


-0.030 


39 


-0.02336 


1.0 


-0.040 


30 


-0.03069 




Table 3.2: CP^ model, charge 1 sector: c and R vs. Ar and At. fo — 1.0 and 
vo = -0.01. 



Ar 


At 


c 


R 


0.100 


0.001 


0.0281 


101. 


0.050 


0.001 


0.0273 


76.6 


0.025 


0.001 


0.0271 


71.0 


0.020 


0.001 


0.0271 


70.3 


0.010 


0.001 


0.0270 


69.4 


0.100 


0.002 


0.0281 


101. 


0.100 


0.004 


0.0281 


101. 


0.100 


0.005 


0.0281 


101. 
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1/[fdot(0,t)]"2 vs ln[f(0,t)] and best fit line to the data 
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-281 1. 69 x+ 10204.2 



-2.5 



-1.5 
ln[f(0,t)] 



-0.5 



Figure 3.4: CP^ model, charge 1 sector: Plot of l/(/)^ vs ln(/) and the best 
fit line to this data. 



Table 3.3: CP^ model, charge 1 sector: c and R vs. Vq. Jq — 1.0. 



vq c R 

-0.01 0.0263 53 

-0.02 0.0485 34 

-0.03 0.0683 25 

-0.05 0.104 17 

-0.06 0.121 15 
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f(0,t) vs. t: Model and Predicted 



o 0.4 




60 

Time in seconds 



120 



Figure 3.5: <CP^ model, charge 1 sector: Predicted course of f{0,t) from 



equation ||3.5|| and actual course of /(0,t) vs. t. 



Table 3.4: CP^ model, charge 1 sector: c and R vs. /q. vq = —0.01. 



/o c R 

1.0 0.0267 62 

2.0 0.0263 108 

3.0 0.0260 150 

4.0 0.0259 190 
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R vs 1/vO and best parabolic fit to the data 
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Figure 3.6: CP^ model, charge 1 sector: R vs 1/vO and the best fit parabola. 
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3.3.2 Characterization of time slices f{r,T) 

Making a closer inspection of the time profiles /(r, T) with T fixed as in Figure 
|3.2| , one may observe that the initial part of the data is close to a hyperbola 
as seen in Figure p7f[ The best hyperbolic fit is determined by a least squares 
method. 

The equation for the hyperbola is 

One would naturally ask about the evolution of the hyperbolic param- 
eters a and h with time, however, neither of these is particularly edifying. A 
simple calculation shows that k should follow f{0,t) closely if b is small, as it 
is. Figure shows an example of a vs. time. Figure shows an example 



of b vs. time, and Figure |3.10| shows an example of k vs. time. In all these 
figures /o = 1.0 and vq = —0.01. 

The evolution of —b/a gives the slope of the asymptotic line to the 



hyperbola, and this evolution is close to linear as seen in Figure p.ll| . One 
might naturally ask how the slope m and intercept bi of the best fit line to the 
evolution of the hyperbolic —b/a change with various initial conditions on /q 



and vq. The variation of m and bi with /o is shown in Table 3.5. It is easy 



to see the slope m varies linearly with l//o, while the intercept bi changes by 



only about 2%. The variation of m and 6j with vo is shown in Table |3.6|. Here 



m and bi vary quadratically with the initial velocity. In Figure |3.12| we show 



the slope m vs. the initial velocity Vq plotted along with the best parabolic 
fit f{x) = -0.11083929x2 + 0.00009768x + 0.00000220 to this data. In Figure 
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|3.13| we show the intercept 6j vs. the initial velocity Vq plotted along with the 
best parabohc fit /(x) = -0.77595779x2 ^ o.00943555x + 0.00006655 to this 
data. 
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Figure 3.7: CP^ model, charge 1 sector: Time slices evolve hyperbolic bump 
at origin. 



Table 3.5: CP^ model, charge 1 sector: Best fit line to evolution of —b/a slope 
m and intercept bi vs. /q. 

/o m bi 

1.0 -0.0000103 -0.000109 

2.0 -0.00000516 -0.000108 

3.0 -0.00000344 -0.000107 

4.0 -0.00000258 -0.000107 
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Hyperbolic parameter a as a function of time 

15 I 1 1 1 1 1 1 1 — 



CO 




20 30 40 50 60 70 80 90 100 110 120 

Time in seconds 

Figure 3.8: CP^ model, charge 1 sector: Plot of hyperbolic parameter a vs 
time, /o = 1.0, vq = —0.01. 



Table 3.6: CP^ model, charge 1 sector: Best fit hne to evolution of —b/a slope 
m and intercept bi vs. vq. 

Vq m bi 

-0.01 -0.0000103 -0.000109 

-0.02 -0.0000436 -0.000426 

-0.03 -0.000100 -0.000916 

-0.05 -0.000281 -0.00235 

-0.06 -0.000402 -0.00329 
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Hyperbolic parameter b as a function of time 
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Figure 3.9: CP^ model, charge 1 sector: Plot of hyperbolic parameter b vs 
time, /o = 1.0, Vq = —0.01. 
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Hyperbolic parameter k as a function of time 
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Figure 3.10: CP^ model, charge 1 sector: Plot of hyperbolic parameter k vs 
time, /o = 1.0, Vq = —0.01. 



Hyperbolic parameter -b/a as a function of time with best fit iine 




Figure 3.11: CP^ model, charge 1 sector: Plot of hyperbolic parameter — 
vs time, /o = 1.0, Vq = —0.01. 
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Hyperbolic parameter -b/a slope vs. initial velocity 
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Figure 3.12: CP^ model, charge 1 sector: Plot of slope m vs initial velocity 
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Hyperbolic parameter -b/a intercept vs. initiai veiocity 
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Figure 3.13: CP^ model, charge 1 sector: Plot of intercept hi vs initial velocity 



Chapter 4 
The CP^ model, charge 2 sector 



It is thought that the charge 2 sector of the CP^ model should behave similarly 
to the 4 + 1 dimensional instanton model. The charge 2 sector consists of all 
rational functions oi z = x + iy of degree 2. A general static solution is of the 
form 

pz + -f 

u = a + — r 

z'' + oz + e 

depending on the five complex parameters a, (3,'j,S,e. To simplify consider 
only solutions of the form 

with 7 real. A straightforward calculation determines the evolution equation 
for the radially symmetric function f{r,t) = 7. It is: 



J J ~^ J, J2 _|_ ^4 ~'~ J2 _|_ ^4 



/2 



(4.1) 



Immediate similarities can be seen with equation p.2|| . The static solutions 
are f{r,T) = c for T fixed and any constant c. We investigate progression 
from /(r, 0) = Cq > towards the singularity at /(r, T) = 0. 
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4.1 Numerics for the CP^ charge 2 sector model 

As with the other models, a finite difference method is used to compute the 
evolution of |[4.1|| numerically. Centered differences are used consistently except 



for 



r 



A naive central differencing scheme here will yield serious instabilities at the 
origin. As before, let 

£/ = r-^oydrf = f" + — . 

r 

This operator has negative real spectrum, hence it is stable. The natural 
differencing scheme for this operator is 

5 



-5 



^ ^ 6Y n{r + 6)-f{r) \ _ f sy f f{r) - f{r - 6) 



2 \ S \ 2 \ S 



This is the differencing scheme used for these terms. 

Now with the differencing explained, we derive f{r,t + At) in exactly 
the same manner as for the 4+1 dimensional model and the charge 1 sector. 
We have an initial guess for f{r,t + At), either given by f{r,t + At) = f{r, t) + 
foAif: with fo the initial velocity given in the problem, or on subsequent time 
steps /(r, t + At) = 2/(r, t) — /(r, t — At). We use this guess to compute 
f{r,t) on the right hand side of [[4.1|| , and then we can solve for a new and 
improved f{r,t + At) on the left hand side of [ [4.1| |. Iterate this procedure to 
get increasingly accurate values of f{r,t). 
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The boundary condition at the origin is found by requiring that /(r, t) 
is an even function, hence 

/(0,t) = ^/(Ar,t)-i/(2Ar,t). 

At the r = R boundary the function should be horizontal so f{R, t) = f{R — 
Ar, t). 

4.2 Predictions 



Equation ||1.4|| gives us the Lagrangian for the general version of this problem. 



We are using 

A 

for our evolution, and via the geodesic approximation we restrict the La- 
grangian integral to this space, to give an effective Lagrangian. The integral of 
the spatial derivatives of u gives a constant, and hence can be ignored. Under 
these assumptions, up to a multiplicative constant, the effective Lagrangian 
becomes 

rdr 







r4 + A2)2 



which integrates to 



A^TT 

L - 



3A 

Since the potential energy is constant, so is the kinetic energy, hence 

A = ky/X. 

Integrating this one obtains 

A = (Cit + C2)^. 
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If A = occurs at time T, we find 

T = 

Cl 

hence we rewrite this as 

X{t)=a{t-Tf. 

Since equation ||4.1|| tends towards the hnear wave equation when r oo, 



the interesting behavior will occur at the origin. This is how we predict that 
/(O, t) will evolve. 

4.3 Results 

The computer model was run under the condition that /(r, 0) = /o with various 
small velocities. The initial velocity is /(r, 0) = vq, other input parameters 
are r^ax = R, Ar and At. 

4.3.1 Evolution of /(O, t) 



A typical evolution of /(O, t) is given in Figure This is best modeled by a 
parabola of the form a{t — T)^ + c, where c is small and grows smaller as the 
grid size decreases. This curve in particular is approximated by 0.0000998(t — 
100.2)2 - 0.00102. This is close to the a{t - Tf predicted. 

We fit /(O, t) to a parabola of the form a{t — T)^ + c for various initial 
conditions. Table ^]l] gives initial conditions /o = /(r, 0) and vq = f{r, 0) and 
the subsequent T and a when Ar = 0.1 and At = 0.005. 

We have 

T ^ — 

l^ol 
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and 



4/0 



Table 4.1: CP^ model, charge 2 sector: Parameters for best fit parabola to 
f{0,t) vs. initial data /o and Vq. 



k 


^^0 


T 


a 


1.0 


-0.01 


204 


0.0000245 


1.0 


-0.02 


102 


0.0000978 


1.0 


-0.03 


68 


0.000220 


1.0 


-0.04 


51 


0.000390 


1.0 


-0.05 


41 


0.000609 


0.5 


-0.01 


104 


0.0000479 


2.0 


-0.01 


404 


0.0000124 


3.0 


-0.01 


604 


0.00000828 




Figure 4.1: CP^ model, charge 2 sector: Evolution of f{0,t) and overlaying 
fit to parabola. 
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4.3.2 Characterization of time slices f{r,T): evolution of a horizon- 
tal ine 

With the evolution of f{0,t) taken care of, we consider the shape of the time 
slices f{r,T) for a given fixed T. This is rather striking, as with the 4+1 
dimensional model, an elliptical bump forms at the origin, as seen in figure 



4.2 



Time slices f(r,T) 




250 



Figure 4.2: CP^ model, charge 2 sector: Time slices of f{r,T) evolve an 
elliptical bump at the origin. 



As with the 4+1 dimensional model, this elliptical bump has equation 
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1. 



(4.2) 



The question naturally arises how the parameters a, b, and k evolve, and they 
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evolve as they did in the 4+1 dimensional model. 



t 



4/0 



k = fo + vot 

Table 0| shows these values as calculated using least squares fitting (to either 



a line or parabola). Since the elliptical fit only works before the right end of the 
ellipse hits the boundary at r = i? the fit sometimes needed to be restricted to 
the portion of the data before this occurred. While the ellipse is small, there 
is a great deal of noise in finding the elliptical parameters, and to get a good 
data fit, this noise must often be removed. Here and is the slope and 
intercept of the line a{t), and likewise rrik and bk are the slope and intercept 
of the line k{t). The parameter c is that in b{t) = ct^ 

A typical evolution of a with /o = 1.0 and vq = —0.01 is in Figure O 
Note the initial noise. 



Figure is a typical evolution for elliptical parameter 6 as a function 
of time with /o = 1.0 and vo = —0.01. 

Figure ^]5| is a typical evolution for elliptical parameter as a function 
of time with /o = 1.0 and vq = —0.01. 
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Tabic 4.2: CP^ model, charge 2 sector, Elliptical parameters vs. initial condi- 
tions /o and Vq 



fo 


^^0 


rria 


ba 


c 


rrik 


h 


0.5 


-0.01 


0.998 


0.050 


0.0000496 


-0.00998 


0.500 


1.0 


-0.01 


0.998 


0.023 


0.0000250 


-0.00999 


1.000 


2.0 


-0.01 


1.000 


-0.138 


0.0000125 


-0.01000 


2.000 


3.0 


-0.01 


0.999 


-0.080 


0.00000836 


-0.01000 


3.000 


4.0 


-0.01 


1.000 


-0.173 


0.00000628 


-0.01000 


4.000 


1.0 


-0.02 


0.995 


0.206 


0.000100 


-0.01994 


1.000 


1.0 


-0.03 


0.989 


0.370 


0.000224 


-0.02983 


0.996 


1.0 


-0.04 


0.983 


0.507 


0.000396 


-0.03965 


0.993 



Elliptical parameter a as a function of time 




Figure 4.3: CP^ model, charge 2 sector: Elliptical parameter a as a function 
of time. 
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Elliptical parameter b as a function of time 




20 40 60 80 100 120 140 160 180 200 

time in seconds 



Figure 4.4: CP^ model, charge 2 sector: Elliptical parameter 6 as a function 
of time. 
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Figure 4.5: CP^ model, charge 2 sector: Elliptical parameter A; as a function 
of time. 



64 

4.3.3 Characterization of time slices /(r, T): evolution of a parabola 

As with the 4 + 1 dimensional model, the evolution of the ellipse suggested the 
curve was trying to obtain the shape of a parabola of the form 



f{r,t) =pr^ + h. (4.3) 

To get the general form of the parabola, we follow the calculation from 
the 4 + 1 dimensional model. From our ellipse equation | |4.2| : 

dy x^lP' 
dx {y — k)a? 

so 

d^y —IP' xb^ dy 

dx {y — k)a'^ {y — kYa^ dx 

At x = 0, y — k = b and this gives 

(fy_-b 

dx 0? 

Recall from the previous section that h = dP and a = t, so this gives 

d-'y 



dx"^ 



The identification of c gives 



So 



dx^ 4/o ■ 

^ 2 dx^ 8/o 



When a run is started with this initial data, /o = vq = —0.02, /o = 

2 

/(0,0) = 1.0 and p = = —0.00005, the time slices of the data have this 



same profile. This is shown in figure [[4.6|| . The parabohc parameter p varies 



by less than 1 part in 10 during the run as seen in figure |[4.7|| . The parabolic 
parameter h evolves close to /(O, t) as seen in figure ||4.8| . 

For a parabola of the form 



f{r,t)=p{ty + h{t), 



we have 



i.e. p{t) is constant, and 



h{t) 



"0 I ^ _ ^ 

1^0 1 



4/o 



Using the identification of p and h in the parabolic form of f{r,t) and 



letting 



we have 



r = t 



2fo 



Substitute this into the partial differential equation |[4.1|| , get a common de- 
nominator, and simplify to obtain 



^0 I r 



mi V 4 



i 



2/o 



6 4 



+ 



^0 4 I ^0 A 



32/3 4 32/3- ■ 2/o- 
The difference between the two sides is 



As with the 4+1 dimensional mode, our concern is with the geodesic approxi- 
mation, and so Vq/I/o) is always chosen to be less than 1/200. The correction 
is then much smaller than the term 



Time Slices f(r,T) 



P 




-3 I 1 1 1 1 1 

50 100 150 200 250 

Time in seconds 



Figure 4.6: CP^ model, charge 2 sector: Time slices of the evolution of 
parabola are parabolas. 



Parabolic parameter p as a function of time 




Figure 4.7: CP^ model, charge 2 sector: Evolution of parabolic parameter 
with time. 




Figure 4.8: CP^ model, charge 2 sector: Evolution of parabolic parameter h 
with time, comparison to f{0,t). 
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Conclusions 

This dissertation investigates the shrinking of sohtons in the 4 + 1 dimensional 
hyperbohc Yang Mills Lagrangian and the charge 1 and charge 2 sectors of the 
model. 

Some theoretical work on the validity of the adiabatic limit for the 
monopole solutions to the Yang-Mills-Higgs theory on Minkowski space is 
presented in 0. 

Of these models, only the CP^ model, charge 1 sector has been inves- 
tigated in any detail in the literature. Many publications 0|, |^, [|], are 
concerned with the translation and scattering of CP^ solitons. The shrinking 
of solitons is investigated in and . The stability of solitons is investigated 
ini. 

In [Q, solitons are found to be numerically unstable. A solution of the 

form 

z 

shrinks spontaneously under their numerical procedure. This does not occur 
in our implementation of any of the three models in this dissertation. The 
static solutions do not evolve in time unless given an initial rate of shrinking. 
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Further, stability and convergence analysis of two of the numerical procedures 
is provided in the Appendix. 

In [0 the time evolution of the shrinking of solitons was studied. They 
arbitrarily cut off the Lagrangian outside of a ball of radius R, to prevent 
logarithmic divergence of the integral for the kinetic energy, and then analyze 
what happens in the R ^ oo limit. In ^ the problem of the logarithmic 
divergence in the kinetic energy integral is solved by investigating the model 
on the sphere 5*^. The radius of the sphere determines a parameter for the 
size analogous to the parameter R for the size of the ball the Lagrangian is 
evaluated on in 

We use a method analogous to that in of cutting off the Lagrangian 
outside of a ball of radius R, and we find an explicit integral for the shrink- 
ing of the soliton, dependent on two parameters: c which is a function of the 
kinetic energy, and R which is the size of the ball on which we evaluate the 
Lagrangian. Once these are specified, this integral gives the theoretical tra- 
jectory of the soliton. The dependence on a cut off R is shown to be a feature 
of the dynamics in the 1+1 dimensional partial differential equation modeled 
in this dissertation, although the dependence of R on the initial conditions is 
unclear. The shrinking predicted by the cut off Lagrangian and the shrinking 
found in the 1+1 dimensional partial differential equation model match. 

In 0, via the cut off Lagrangian, an ODE for the evolution is found, 
and approximate solutions to this ODE were evaluated. A 1+1 dimensional 
partial differential equation was also solved numerically in 0] for the evolution 
of the shrinking of solitons, although the initial condition is not a static soliton. 



They find the shrinking rate is given by a power law corrected by a small 
logarithmic term, with the corrections vanishing as — > oo. In contrast, we 
find an explicit integral for the shrinking of a soliton where the Lagrangian is 
cut off outside of the ball of radius R and that the corrections between this 
and the R ^ oo limit are necessary. 

Similar analysis of the predicted shrinking is given in the CP^ model, 
charge 2 sector and the 4+1 dimensional model. These two models have no 
divergences. A prediction is made by evaluating the effective Lagrangian, and 
the evolution of the 1+1 dimensional partial differential equation is shown to 
be close to that predicted. 

In addition to this, in all of the models, the shape of a time slice /(r, T), 
with T fixed, is characterized, by elliptical bumps at the origin in the 4+1 
dimensional model and the CP^ model, charge 2 sector or by hyperbolic bumps 
at the origin in the CP^ model, charge 1 sector. 

In the 4+1 dimensional model and the CP^ model, a better approxima- 
tion than the geodesic approximation for solutions to the partial differential 
equations are found in the evolution of the parabolas. 



Appendices 
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Appendix A 



Stability and convergence of the 4+1 
dimensional model 



In this chapter we analyze the stabihty of the equation 



r 



2/2 2{ff 8rf 



(A.1) 



f + r"^ f + t"^ f + t' 



.2 ' 



and the associated differencing scheme used in finding numerical solutions. 
Here / = f{r,t), and r is a radial variable, hence r > 0. It will be shown that 
the stability behavior of the differential equation is qualitatively the same as 
that of the difference equation, and that the difference equation converges to 
the differential equation as Ar — > and At — > 0. 

The usual Von Neumann stability analysis has one substitute f{r,t) — 
fo{r,t) + £e"^*e*'*'', where fo{r,t) is presumed to solve the equation exactly, 
linearize in e and then solve the resulting equation for a; in terms of k,. In this 
case, if cu has a negative imaginary part, then it has a growing mode. The 
equation is stable if it has no growing modes. This analysis is only valid when 

K ^ z and -. 



length scale 



r 



74 



75 

A.l Continuum stability: simplified model 

The first thing to address is the stabihty of the linear part of the partial 
differential equation [|A.1||. This is: 



r 
Set 

and plug into |[A.2|| to obtain the equation: 

U'' = 

r 

Solving this for u yields: 



/ = f" + — . (A.2) 



^ / 2 ^^'^ 
UJ = ±\ K"^ 



So, by this analysis there is always a solution for uj with a negative 
imaginary part, and therefore this equation has a growing mode and is not 
stable. 

However, if we consider abstractly the equation 

f = f" + — = ^f 
r 

where Lf is a linear operator; it is easy to see that 
and hence that 
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where B = r^^^ and 



A^A is hermitian, and so it has real spectrum. Because it is essentially a 
square, it has positive real spectrum. Consequently spec(£;) = spec (—^4^ A) is 
real and negative. So the solutions to the equation 



consists of sines and cosines in the time variable multiplying the eigenfunctions 
of £, and hence has no growing mode and is strictly stable. 

So what gives? 

The problem stems from the Von Neuman stability analysis which sets 
f{r,t) = e*''^e*'^* instead of using the more general stability analysis where 
f{r,t) — g{r)e^'^*. The Von Neumann form is used because it results in al- 
gebraic equations instead of differential equations, which are much easier to 
solve. The Von Neumann analysis yielded 



Since in this analysis 1/r, approximate using Taylor's theorem to obtain 

, hi 

UJ ~ ±K ^ — . 




2r 



Since velocity v is 



du 



OK 



then r To ^ so 
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5 



What happened was that we let g{r) = e*'"^' and so the r behavior had no 
choice but to come out in the e*"^* portion of the equation. 

If, however, we forcibly put the factor of r~^/^ into the Von Neumann 
stability analysis, via 

/(r,t) = r-5/2e*'^'-e'^* 

plugging / into ||A.1|| one obtains the equation: 

15 



2 2 I 

- = - +4r2 



which has solutions. 



1 VT5 + 472k2 
cu = ± 

2 r 



So now, uj has only real roots, because the "correct" form for g{r) was used. 
A. 2 The whole equation 

Similarly, any analysis of the general equation ||A.1|| will have some "growing 



modes" analogous to the factor of r~^/^ found in the linear equation that do 
not affect the general stability. One can put this factor in explicitly to remove 
these modes. 

Set 

/(r,t) = /o(r,t)+5r™e*'^V^* 
and linearize ||A.1|| in e to obtain the equation: 
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2 4iu;/o 
K — ^ + --^-f 



15 



_8r^ 10(/^ + 2r) 

/o + r2 ' 4r2 ' (/o + r2)2 r(/o + r^) 

| 4i^(/^ + 2r) ^ 2(/^ + 2r)^ ^ 2fl 



fo + r-" 



Solutions to this equation are found to be 



UJ 



= ^(-4ir/o 



± 



8r Vo + 4«V' 



2„6 



1 1/2n 

-SO/or^ - 65r^ 



Now address the realm where fo and /q are small, and k is large. Use 
a first order Taylor approximation to the square root to obtain, 



-2i/o 



(/o + 



1 + 



2z (/^ + 2r) 5[,/o /o^ - /^^ 



72r/^ + 50/o + 65r2^^/2\ ^/ 1 



^4 J, 2 



(A.3) 



There arc clearly some negative imaginary parts of a;, but since k S> - they 

r 

are all bounded. 
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A. 3 Discretization scheme 

After much preliminary work on stability with a standard finite differencing 
scheme, it was discovered that the main stability problems were generated by 
the linear part of the equation, namely 



r 

as r — > 0. Using a naive centered difference scheme for this part of the equation, 
one has irreconcilable problems with an instability the origin which is not 
present in the continuum equation. The "general rule of thumb" to use is 
when you have a problem, discretize the problem in the natural way for its 
differential operator. Proceeding along these lines: 



2)/ = r + — 

r 

and it was shown in section 2.1 that 2) = r^^drV^dr- Discretize with q indexing 
the space variable r so qAr = r and n indexing the time variable t so nAt = t. 
The discretization for 2)/ is 



{qAr) 



q + -]Ar 



f{{q + 1) Ar, nAt) - /(gAr, nAt) 



Ar 



Ar 



/(gAr, nAt) - /((g - 1) Ar, nAt) 
Ar 



Ar 



This differencing scheme removed the problems at the origin completely. 
Other than this, centered differences are used to discretize the equation ||A.1[ . 
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A. 4 Outline of the stability analysis 

The stability analysis for this discretization scheme will be approached in the 
following way. 

First, analyze the stability as r — > by linearizing the difference scheme 
as a matrix XL and finding its eigenvalues and eigenvectors. 

Second, use Von Neumann analysis on the difference scheme in the 
realm where r is bounded away from zero and /(r, t) = /o(r, t) + £e*'"^'^''e*'^"'^* 
where fo{r,t) = c, a constant; fo{r,t) = c does solve the equation exactly. 
Find cu in terms of k and show that the growing modes are all bounded. 

Third, analyze what happens when fg is not presumed to be constant, 
but the derivatives are presumed to be close to zero. The equation for cu in 
terms of n under this new circumstance gains additional dependencies on k 
and uj. Given k, one can show that inverting this map on cu yields contraction 
maps for cu in balls of radius \k\/2 about cu — ±k. So by the Contraction 
Mapping Principle, there are solutions for cu in the balls about ±k. From the 
form of the equation it can be seen that there are two solutions for uj, therefore 
these are they. Now let cuq indicate the solution for cu with /o presumed to 
be constant, and cui indicate the solution for cu with no such assumption, one 
can use the conclusions of the Contraction Mapping Principle to show that 
the difference \u!i — cuo] is bounded. Hence the negative imaginary part of cui, 
i.e. the growing mode, is bounded. 



A. 5 Stability near zero 

The Von Neumann stability analysis is not valid for r — 0. To analyze the 



stability here, let / = fo + e6f and hnearize the equation ||A. 1|| in e. Represent 
this as 

6f = L,i6f) + L2i6'f), 

where £i and £2 are linear operators, and £2 is close to a multiple of the 
identity matrix, and so one can replace it with cl. Then discretize in space, and 
find the eigenvalues and eigenvectors of £1, and relate these to the eigenvalues 
and eigenvectors of 



( I\ 


' 5f - 






cl J 









(A.4) 



This equation is 



Mv = v. 



if A is an eigenvalue then one has 

V = \v. 

We know the solutions to this, for the i*"^ component of v we get 

Vi = Ae^\ 

Hence Vi does not grow with time if A does not have a positive real component. 

Now, to relate the eigenvalues and eigenvectors of £1 to the eigenvalues 
of the matrix in ||A.4||, presume that the vector 



a 
b 
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is an eigenvector of the matrix in equation ||A.4|| with eigenvalue A. Therefore 



( ° 




a 


= A 


a 


V ^1 


i) 


b 


b 



The top part of this equation states that b = Xa. Now one can use this in the 
bottom part of the equation in the following manner: 

£ia + cb = Xb 

£jia + cXa = }?a 

In particular one sees that a must be an eigenvector of £i, and that if the 
eigenvalues of £1 are equal to a then one can find the eigenvalues of the whole 
matrix via the equation 

aa + cXa = X^a. 

The solutions are 



c± Vc2 + 4a 



So if, as we are about to find, a < 0, we have two cases: 



c < : 3?(A) < 
: 3?(A) < c 



So any solution of [A.4| can be written as a linear combination of solutions 
with at worst bounded growth rates. 

Now let's turn to showing that a < if a is an eigenvalue of £1. Allow 
6f{q, n) to represent 6f{qAr, nAt) and likewise with fo{q, n) = fo{qAr, nAt), 
compute the following linearized equation, which is discretized in space only: 



83 



-2 



' 2g (5/ g + 1, n) - 5/ g - 1, n) - — . . . ..^ ^ 



+(gAr)- 



g + - ) Ar 



/o(g,n) + (gAr)2 

Y5/(g + l,n)-5/(g,n; 



Ar 



g Ar 

^ 2 ' 



Ar 



Ar 



+4 



{qArfSf{q,n) 
(/o(g,n) + (gAr)2f 

fo{q,n)6f{q,n) 
Mq,n) + (gAr) 



(/o(g,n))^^/(g,n) 
(/o(g,n) + (gAr)2f 



One is concerned about behavior as r ^ 0, and so the functions 
/o(gAr, nAt) and /o(gAr, nAt) can be approximated by /o(0, nAt) and /o(0, nAt) 
which shall be called /o and fg. Likewise, /o + (gAr)^ will be replaced with 
/o, and this yields the following equation: 



^'fiq,n) = 



2q{5f{q + l,n)-5f{q-l,n))- 



Aq\Arf6f{q,n) \ 
fo 



fo 



+(gAr)" 



g + - ) Ar 



5f{q+l,n)-5f{q,n) 



Ar 



q Ar 

^ 2 ' 



{qArf5f{q,n] 



Sf{Q,n) -Sf{q-l,n] 
Ar 

{forSf{q,n) 



Ar 



(A.5) 



ifoY 



-2 



ifoY 
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+4 



\ fo 



n] 



Since this analysis is for r — > 

J 

fo 



^2 = 



Cl 



is a good approximation. 

To find £i use ||A.5|| and the quadratic fit boundary condition at the 
origin, i. e. 

6f{0,t) = pf{Ar,t)-^6f{2Ar,t). 

If one lets 

and fo{k,t) = fo{kAr,t), then one obtains the following tridiagonal matrix: 



+ (Ar) 



5 /^\ 5 



fo 
fo 



ai,2 



(Ar) 



^ + k'\Ar)-' 
Jo 



2 
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0'k,k+l — 



-Ak 



k+-\ -[k- 



- 2 



/o 



/o 



^k-^Ar) 



One can use Maple to compute the eigenvalues and eigenvectors of this 
matrix while changing the size of the matrix n, and the values of /o, /o, and 
Ar. 

Since /o is always much less than /o, the contributions from terms with 
/o are negligible. Further, reparametrizing by multiplying a factor p times /o, 
Ar and times /o yields eigenvalues times the originals. In fact, the 
largest of the terms in any of these matrix elements should be that multiplied 
by Ar~^, hence it is to be expected that rescahng Ar itself by a factor p should 
result in the eigenvalues changing by approximately 1/p^. 

For example, when n = 5, /o = 1, /o = —0.01, Ar = 0.01 the matrix 



is: 



-75828. 75828. 
2381.1 -32891 











30510. 
4030.7 -25633. 








5145.1 







21602. 

-23149. 18004. 

5925.0 -22010. 



It has eigenvalues with multiplicity m and eigenvectors given by: 



Eigenvalue m 

-79876. 1 

-42162. 1 

-31531. 1 

-18654. 1 

-7316.6 1 



Eigenvector 

8.8171, -.47042, .036302, -.0034019, .00034922.] 
-2.0274, -.90138, .43183, -.16212, .047685] 
6.3939, 3.7401, -.33103, -.60680, .37688] 
1.7172, 1.2964, .47128, -.08927, -.15705] 
-1.6083,-1.4544,-1.0938, -.65654, -.26464] 



Now rescale by p = 5 and use parameters /o = 25, /o = —0.05, Ar 



0.05, n = 5 to get the matrix: 



-3033.1 
95.243 









3033.1 
-1315.6 
161.23 









1220.4 
-1025.3 
205.80 








864.08 
-925.96 
237.00 








720.17 
-880.40 



This has eigenvalues with multiplicity m and eigenvectors: 



Eigenvalue 
-3195.9 
-1686.6 
-1261.1 
-745.84 



m 

1 
1 
1 
1 



-292.56 1 



8.8181, -.47043 
-2.0282, -.90108 
6.3888, 3.7349 
1.7173, 1.2958 
-1.6088,-1.4547 



Eigenvector 

.036302, - .0034009, .00034900] 
.43197, -.16222, .047711] 

-.33118, -.60692, .37736] 
.47102, -.08941, -.15695] 

-1.0937, -.65607, -.26455] 



These are 1/25*^ of those of the previous case, as expected. 

Now, vary /o, using parameters /o = 1, /o = —0.1, Ar = 0.01, n 
and get the matrix: 



-75828. 
2381.1 






75828. 

-32891. 
4030.7 







30510. 
-25633. 
5145.1 







21602. 
-23149. 
5925.0 







18004. 
-22010. 



With eigenvalues: 

Eigenvalue m 

-79876. 1 

-42162. 1 

-31531. 1 

-18654. 1 

-7316.6 1 



Eigenvector 

8.8171, -.47042, .036302, -.0034019, .00034922] 
-2.0274, -.90138, .43183, -.16212, .047685] 
6.3939, 3.7401, -.33103, -.60680, .37688] 
1.7172, 1.2964, .47128, -.08927, -.15705] 
-1.6083,-1.4544,-1.0938, -.65654, -.26464] 



These are unchanged to this level of precision from those found with /q 
—0.01. The contribution of the terms with /o is negligible, as expected. 
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Now, rescale Ar = 0.1 while leaving /q = 1 and /o = —0.01 and get 
the matrix and eigenvalues: 



-753.00 753.00 

31.731 -328.91 

52.187 







297.18 

-256.33 204.14 

67.291 -231.49 164.20 

79.050 -220.10 



Eigenvalue m 

-806.30 1 

-436.27 1 

-317.81 1 

-175.73 1 

-53.872 1 



Eigenvector 

8.1450, -.57701, .057337, -.0069800, .00094304] 
-2.3858,-1.0045, .61732, -.28715, .10502] 
4.0282, 2.3296, -.34264, -.49247, .39909] 
1.3007, .99754, .37537, -.10673, -.19021] 
1.2249, 1.1373, .92190, .62377, .29669] 

So rescahng Ar by a factor of 10 rescales the eigenvalues by approximately a 
factor of 1/100^^, with the scaling best closest to the origin. This is consistent 
with the original equations. 

Lastly, check to see what effect rescaling n or the size of the matrix has. 
Ideally since we want to explore behavior near r = we would want to take 
n — > oo as Ar — > and nAr — > 0, but this isn't possible, so instead, simply 
double the size of n to n = 10. The matrix and eigenvalues are: 



-75828. 75828. 
2381.1-32891.30510. 
4030.7-25633.21602 












5145.1-23149.18004. 

-5925.0 -22010. 16085. 



6496.2 -21394. 14897. 



6931.7-21023.14091 









7273.9 -20783. 13509. 

7550.2 -20618. 13068. 
7777.8 -20501. 



Eigenvalue m Eigenvector 



-79898. 


1 


[8.8170, 


-.47043, 


.036302, 


-.0034177, 






.00035827, 


-.000049446, 


.35708 xlO-^ 


-.44969 xlO-^ 


-42854. 


1 


[2.0301, 


.88204, 


-.44662, 


.19165, 






-.082308, 


.036070, 


-.016080, 


.0071721, 






-.0030613, 


.0010666] 










-.6539 xlO-^ 


-.92087x10-6] 






-39076. 


1 


[7.1614, 


3.4681, 


-1.2625, 


.13923, 






.23670, 


-.30248, 


.25528, 


-.17811, 






.10369, 


-.043440] 






-34796. 


1 


[9.2387, 


4.9953, 


-1.0334, 


-.49271, 






.61276, 


-.30607, 


.00732, 


.14346, 






-.15312, 


.08324] 






-29191. 


1 


[-19.285, 


-11.852, 


.068578, 


2.1976, 






-.75346, 


-.47022, 


.57537, 


-.10236, 






-.24492, 


.21945] 






-22916. 


1 


[9.1374, 


6.3726, 


1.3700, 


-1.0152, 






-.40642, 


.39472, 


.13659, 


-.21252, 






-.04067, 


.13029] 






-16554. 


1 


[-3.3235, 


-2.5963, 


-1.1309, 


.00853, 






.32669, 


.10839, 


-.10703, 


-.08724, 






.03047, 


.059846] 






-10673. 


1 


[8.5362, 


7.3303, 


4.6724, 


1.8693, 






-.040111, 


-.71838, 


-.50033, 


-.014049, 






.25853, 


.20488] 






-5754. 


1 


[1.7509, 


1.6167, 


1.3015, 


.89606, 






.49400, 


.16895, 


-.03811, 


-.12426, 






-.11777, 


-.062081] 






-2158. 


1 


[.57440, 


.55752, 


.51684, 


.45759, 






.38593, 


.30774, 


.22905, 


.15526, 






.090760, 


.038476] 







As before, all eigenvalues are negative. The modes localized near zero don't 



change. The eigenvalues increase towards zero as n increases because there 
are growing modes of this equation away from r — 0. 

One can conclude that the eigenvalues of this matrix under reasonable 



initial conditions will always all be negative, as required. This differencing 
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scheme does not have any instabihties generated at the origin. 



A. 6 Stability away from zero: simplified model 



In this section the stabihty of difference equation derived from ||A.1|] with 



the differencing scheme outhned in section 2.1 under the special condition 
that f{r,t) = fo{r,t) -)_ 5e*'"?'^''e*'^"^* where fo{r,t) = c, a constant, will be 
addressed. 

An astute reader will note that the factor of r~^/^ is missing in this 
analysis. Its explicit inclusion makes the resultant expression for u far uglier 
than even it is now. A part of the growing mode is analogous to the r~^/^ 
term, but not exactly the same. 

Plug into ||A.1|| and linearize in e to obtain the following equation: 



+ (e^'^^'' - 1) - - (l-e-«^^) 



(At)2 g5(Ar)2 

/o + (gAr)2 • 

Let X = e'^^^ and 



(A.6) 



"^^ /o + (gAr)2 ^ g5(Ar)2 

Hence J represents the right hand side of the equation. One can reduce equa- 
tion ig to 

a;2-(2 + JAt2)a; + l = 
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and so it has two solutions if we try to solve for x = e*'^*. Hence it has two 
(logarithmic) solutions for cu. These are obtained via the quadratic formula 



so 



-^ln I 1 + ^— ±VJAtWl + ^p 



u — 



At 

Clearly the biggest concern is determining exactly how big JAt'^ is. 
Let kAt = 6, then reduce exponentials to sines and cosines appropriately, and 
expand the factors of g±l/2, and simplify using such information as g = r/Ar, 
to obtain 

— 8mr /sin^\ r,/2cos^ — 2\ 5m /sin ^ 
J = ^— ^ ^ + :r. + 



fo + r^ \ J \ J r \ e 

5K^Ar^ /2cos^^ — 2\ 5m Ar^ /sini 
+ . . + 



2r2 V J 2r^ \ e 

^bn'^Ar^ /2cos^-2\ ^ mAr^ fsinO 



Analyze the behavior as Ar — > and At 0, under the assumption that r > 
is fixed and that At -C Ar. One also knows that kAt < tt. Consequently, all 
terms go to zero as Ar — > except the following: 

-SiKr fsm9\ ^ pcos^ ~ + ^^'^ /^sin^\ ^ bn^Ar^ f2cos9 - 2 



fo + r^ \ J V J r \ 9 J 2r^ \ 9'^ 

The trigonometric part is strictly bounded, and the simple condition that 
At <S cAr with c chosen sufficiently small, implies that KAt^ or [nAtY <^ 
(cnY can be made as small as one likes. 
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In particular, one may use the Taylor approximation to the square root. 
Then one may then use the Taylor approximation to the logarithm, and with 
some cancellation between the first and second term, one obtains: 



-i In 



At 



±VJAt + O ({^/JAt) 



At 



Recall, as Ar — > all of the terms in J go to except for 



-8mr / sin 9 



/o + r2 V e 

So one has 



2cose-2\ 5m /sin ^\ Sre^Ar^ /2cos^ - 2 



—Sinr f sm9\ „/2cos6' — 2\ din /sin 9 

+ TT. + 



/o + r^ 



e 



r \ e 



+^r^ + O Ar 



2r2 



^2 



Since k is large while Ar is small, the largest term here should be the 
unmultiplied by Ar, and so approximating the square root this can be 
modified to: 



IK 



2-2cos^ 



^2 



1 + 



Air 



sva.6 



«(/o + r2) V ^ / V2-2cos^ 
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5i ( sin d 



2Kr \ 



9^ 

2 - 2cose 



4^2 



+ 



Therefore the leading order term of y/j is of size 0{k,). RecaUing nAr < tt 
one can make the correction O as small as one hkes by requiring 

that At < cAr^/^ and choosing the factor c appropriately. And given this 
choice of At, then one has that 



± 



K 



2 - 2cos^ 
Air ( sin Q 



^2 



5i / sin 6' 



terms that go to zero 



(fo + r^)\ J\2-2cos9J 2r \ 9 J\2-2cos9 

bounded terms 



Furthermore it is clear that if one permitted r — > that the dominant piece 
would be 



2 — 2 cos 9 hi f sin 9 



9^ 



2r\ 9 



9' 



2 -2cos9 



and as ^ ^ this goes to 



hi 
2^ 



or exactly what we obtained in appendix A.l that was so often "corrected" by 
the addition of a factor of r~^/^. 



A.7 The complications 
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Now we must deal with the reahty that we have thus far assumed /o = c, and 
in general this is not so. We need to know that when /o is a function with 
small bounded derivatives, that the growing modes remain close to those in 
the previous analysis. 

In the previous section we let f{r,t) = f{qAr,nAt) = fo{r,t) + 
^^iKqAr^tumAt ^ndcr the assumption that /o = c. Call the solution for lu in 
terms of k in this analysis cuq- The equation solved to do this was 



= + 



(At)2 q^{Ary 

^q^Q-iKAr _ gjKAr^ 

/o + (gAr)2 

Characterize this equation as 

g{u) =old(K). 



If one does not assume that /o = c then one gets the equation 

5 / 1 \ 5 

= + 



^g^g— mAr giftAr^ ^jf / ^inAr g— i/tAr 



/o + (gAr)2 /o + r2 V 2Ar 



(/o + r2)2 (/o + r2)2 (/o + r2)2 
4/o /e''^^* - e-^'^^*' 



/o + r2 V 2 Ai 



(A.7) 
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This equation is quadratic in e*'^^*; hence it has two solutions for e*'^^*, and 
hence two solutions for uj for each branch of the logarithm. We wish to get a 
handle on these two solutions for uj. Equation |[A.7|| can be characterized as 



g{ijj) = old(fi:) + newi(fi:) + new2(ci;) 



(A.8) 



with 



newil^Kj 



4/^ fe 



— e 



2Ar 
2/0^ 



(/o + r 



2\2 



and 



new2(u;) 



(/o + r2)2 (/o + r2)2 
4/0 fe 



(A.9) 



(A.IO) 



/o + r2 V 2 At 
Here /o is substituted for the finite difference 

/o(gAr, {n + 1) At) - /o(gAr, (n - 1) At) 
2At 

and likewise with /q and the centered difference equivalent in the r variable. 
Now note that 

, , 2coswAt-2 ,^ ^ , 

9i^) = — . (A.ll) 



and 



and 



At2 



ujo = g ^(old(K)) 



uJi = g (old(/t) + newi(/t) + new2(co'i)). 



Further if one lets 



ujg = g ^(old(K) + s(newi(K) + new2(co's))) 



95 



then 



\UJI - LOq\ 



ds 



ds 



< \[{9 )'](newi(K) + new2(u;J)| (is 

^0 

< max (K^f^"*^)')!) max (|newi(ft) + new2(tt's)|) (A. 12) 



So finding a bound on uj in terms of k will allow one to find a bound for 
ItUi — cuol in terms of n. 

Now we will show that the map 

T(uj) = 5'~^(old(fi;) + newi(K) + new2(u;)) 



is a contraction mapping on the ball — /t| < \k,\/2 = B{k, k/2). Using ||A.ll 
one finds that 

-2sin(a;At) 



At 



We've already stated one needs to choose At < cAr^/^ so specifically, one can 
make At < Ar/150. Then since |/tAr| < vr and Ic^l < 3|k|/2 this implies that 
IcijAtl < tt/IOO and one has 

sin(ti;At) 



\g'{uj)\ = \2uj\ 



> 1.98L' . 



Now if uji and 002 are in B{k, k/2) then assign yi and 2/2 as follows and 
find that 

Ui = old(fi;) + newi(ft) + new2(ci;i) 
y2 = old(K) + newi(K) + new2(ci;2) 



Z/2-Z/1 = new2(a;2) - new2(a;i). 
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Letting, 



calculate 



y{s) =yi + s(y2-yi), 



|T(a;2)-T(a;i)| 

\9-\y2) - 9-\yi)\ 



{9 ) {y{s)){y2 - yi)ds 



< / max\{g ^)'{y{s))\\{y2-yi)\ds 
Jo 

< max 1(2/2 (A.13) 



Simple calculus yields: 



9'i^) 



By choice of yi and y2 and wi,W2 G B{k., k./2), 

1 



max ) 



< 



1.98|u;| - 1.98\k\ 
Since |cc;i|, |tc;2| < 3|fi;|/2 and |ci;iAt|, |tc;2At| < tt/IOO, we have 

4/0 /i sin(c<j2Ai) isin(a;iAt) 



(A.14) 



|new2(a;2) - new2(cui)| 



fo + r' 
4/0 



At 



At 



fo + r' 



f sin(a;2At) sin(a;iAt) 



V At 



At 



< 




cos{sAt)ds 



uj2-u;i . 



(A.15) 
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Now plug [pTTl and [pTTSf into [l/TT^ , to obtain: 



\T{u^)-T{uj,)\ = \g-\y^)-g-\y,)\< 



1.98k| 



4/0 



(/o + r2 



- ^1 



The quantity 



/o 



defines the time scale for the problem. In order for this entire analysis to make 
sense one expects 



K 3> time scale and k ^ length scale. 

Therefore, this is sufficient to show that T{w) is a contraction map from 
B{k,,k,/2) to itself. So by the Contraction Mapping Principle, one concludes 
that there exists a fixed point of T in B{k,,k,/2). Clearly this argument 
holds just as well for — T and B{—k,,k,/2), so there is a fixed point for — T 
in B{—K, '*/2). Since the left hand side of |[A.7|| is quadratic in e*"^^*, there are 
two solutions for u and these are they. 

Lastly, we finish estimating in equation [|A.12|| . We already have a 
perfectly good estimate for (g^^)' from |[A.14|| . We merely need to estimate 



|newi(/t) + new2(t^)| . 



Replacing the exponentials with their trigonometric forms in | A.9 | and [ A. 10 
to get: 



newil^Kj 



4/^ /2^sin(«:Ar)\ , 2 (/^f 



2Ar 



+ 



(/o + r 



2^2 



2/0^ 
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and 



new2((x; 1 



4/0 



2i sin u; At 



/o + rH 2 At 
Recall also that fi;;»l/rsoKr:^l. Now if 9 is real, 



(A.16) 



Isin^l < e, 

and if 6 is complex with \6\ sufficiently small, as it would be if = uAt, then 

Isin^l < l.Oll^l, 



one has 

|newi(K) + new2(co')| < 



4/^ 



2(/^) 



ifo + r 



2^2 



+ 





2/0^ 




4/0 


+ 


(/o + 


+ 


fo + r' 


By my previous arguments about 


K, and the 


len{ 


5th and 



(/o + 
1.01u;| 



central terms are much much smaller than the others. Plugging |ti;| < 3|k|/2 

ne c< 

4/^ 



and this and [|A. 14|] into ||A.12|| one concludes that under these assumptions 

2(/^ 



l.QSkl 



fo + r^ 



+ 



2/0^ 



(/o + r 



2^2 



1^ + 



4/0 



3.03k 



(/o + 



/o + r^ 



Hence, \uji — ujq\ is bounded, which implies in turn that the imaginary part 
of uji is bounded since the imaginary part of c^o was. Once again the growing 
modes are bounded. 



Further, if we compare with equation [|A.3|| , we see analogues between 
the terms in these two equations. 



1.98 



4/^ 



2^(/o + 2r) 
fo + r' 
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and 

3.03 4/0 -2z/o 
1.98 /o + r2 ^ /o + r2- 

A. 8 Convergence estimates 

This section contains an analysis of the convergence of the differencing scheme 
for the equation 

/ = r + ^ + ^-^-^. (A.17) 

^ ^ I ^2 ^ I ^2 I 1^2 

We will show that as Ar — * and At — > that the solution found for f{r,t) 
converges to an actual solution of the partial differential equation. 

First consider the scheme to forward integrate the equation and the 
error in it. Substitute the appropriate differences in ||A.17|| then solve for 
/(r, t + At) in the difference for / i.e. solve: 

/,(r,( + A() = 2Ur,t)-Ur,t-At) 



r fc + r"^ fc + r'^ fc + r"^ 
Here fc is used instead of / to indicate that this is the numerically calculated 
value rather than the one that is a solution to the PDE. Also assume the deriva- 
tives are represented by the appropriate differences. Approximate /c(r, t + At) 
for use in by using first fc{r,t + At) = 2fc{r,t) — fc{r,t — At), and then 
one iterates 6 times, finding a new fc{r, t + At) from solving the equation with 
the previous value. Using greater numbers of iterations on fdfjt + At) does 
not change the answer to the numerical precision on the computer, so this will 
not be considered as a source of error. 
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The error comes from the discretizations. Let At^Ci be the error in 
the difference used for /, and C2 be the accumulated error in the differences 
for the left hand side of |[A.17|| , then one can see that the error is of the form 



\f{r,t + At)-f,{r,t + At)\ = | At^Ci + At^Csl < At^ (|Ci| + IC2I) . 

The error in going 2 steps would be 

\f{r,t + 2At)-Ur,t + 2At)\ < 2At'' {\C^\ + IC^]) 

+e^^* [Af {\C,\ + \C,\)] . 

With the term involving e^^^ coming from the stability analysis, with W > 
\^{iuj)\. One knows from the stability analysis that the negative imaginary 
part of LO is bounded, hence W is. The error in going n steps would therefore 
be: 

\f{r,t + nAt)- f,{r,t + nAt)\ < nAt^ (|Ci| + IC2I) 

n 

+ J2^wjAt [At' {\C^\ + \C,\)] . 
i=i 

And one also knows that 



y^w^At^ / e'^^^'dj = ^———^ (A.18) 

^ - A WAt ^ ' 

Since W was strictly bounded, one represents the above as Now see that 
the error in going n steps is 

|/(r,t + nAt) -/,(r,t + nAt)| < nAt^ (|Ci| + IC2I) 

+AtC3(|Ci| + |C2|). 
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Since C3 is bounded, one merely needs to show that Ci and C2 are bounded 
to know that as At —>■ clearly this error goes to zero. 

Now we use Taylor's Theorem to find the errors in using the differences, 
i. e. to find Ci and C2. 

We use a generic function h{x) and approximate to fourth order: 



Here Ci & [x, x + 5], (2 ^ — S, x]. Use a third order approximation, subtract- 
ing, plus a little algebra, to obtain: 



where ^1 G [x - S,x + 5] and replaced h"'{Ci) + h"'{C2) = 2/i"'(6)- Likewise, 
using the fourth order approximation, add the equations, and do some algebra 



h{x + 6) = h{x) + 5h'{x) + ^\"{x) + ^\"\x) + ^/^('HCi) 

Z D z4 

h{x-5) = h{x)-6h'{x) + ^\"{x)-^\"'{x) + ^h^'\C2). 

2 o 24: 



h{x + 6) - h{x - 6) 52 
26 6" 



h'"i^,) = h'{x). 



to obtain: 



h{x + 5) + h{x -6)- 2h{x) _ 6^ 

52 



h^'\^2) = h"{x). 



Applying this to equation |[A.17|| , one has 



/(r, t + At) + f{r, t - At) - 2/(r, t) At' d^f 
At2 12 dfi 



(r, r]i) 



2p{r,t) _ ifir,t + At)-f{r,t-At)f 



f{r,t)+r' 2At2(/(r,t) +r2) 
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2(f(r,t))^ 



{f{r + Ar,t)-fir-Ar,t)f 
2Ar2(/(r,t) +r2) 

/(r + Ar, t)- f{r- Ar, t) \ / Ar^ 9^ 



Ar (/(r,t) +r2 
2 /Ar^aV 



6 



8rr(r,t) 



f{r,t)+r^ V 6 ar3^^''^^ 
4 (/(r + Ar, t) - /(r - Ar, i)) 



Ar (/(r, +r2) 
8r /Ar^aV 



/(r, + r^ \ 6 9r^ 
The first of tliese equations says tliat if 



(r, 7^1 ) 



is bounded then the error Ci goes to zero as At'^ — > 0. 

If /(r, t) is chosen in C(4, 4), the space of functions with four bounded 
space/time derivatives, all of these differences have bounded errors. 

Now attending to the discretization for 



f"{r,t) + 



5f'{r,t) 



r H — Ar 
2 



1 ^ 

r Ar 

2 



f{{r + Ar,t)-f{r,t) 
Ar 

f{r,t)-f{r-Ar,t) 
Ar 



(A.19) 



(r^Ar), 



one needs to find an approximation of the error on this. Start with h{r) — 
T^f'ir). Consequently: 

h\r) = 5rV'(r) + r'f"{r) = (f"{r) + 



h"{r) = 20rV'(r) + 10r7"(r) + r7'"(r) 

h"'{r) = 60rV(r) + 60rV"(r) + 15rV"(r) + r^'f^^^r). 
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By taking a Taylor series for h{r), obtain: 



Ar 



Ar 



Ar^ 



Ar^ 



h[r + —j = h{r) + —h'{r) + ^h"{r) + ^h"'{C,) 

h{r-—] = h{r)-—h\r) + ^h"{r)-^h'"{C2). 
Subtracting and doing some algebra as before yields: 



h'{r) 



Ar 24 
Now start substituting in definitions for h{r), h'{r), and let 

Ar Ar 



P3 e 



r , r H 

2 ' 2 



etc. to obtain: 



r' f"{r) + 



5f(r) 



(r + ^)V'(r + t)-(--f)V'(r-f) 



Ar 



Ar^ 



(60rV'(p3) + 60rV"(p3) + 15rV'"(p3) + rV^'nP3)) • (A.20) 



A centered difference approximation gives: 



f r + 



Ar\ f{r + Ar) - f{r) Ar' 



Ar 



24 



r(6), (A.21) 



and similarly 



fir-- 



Ar\ f{r) — f{r — Ar) Ar'' 



Ar 



24 



r(6) (A.22) 



So substituting ||A.21|| and |[A.22|| into ||A.20|| , one has the difference equation 



and approximation for [|A.19|| . Letting G [r, r + Ar] and ps G [r — Ar, r] 
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r + 



Ar 



1 5 



Ar 



f{r + Ar,t)-f{r,t) 
Ar 

f{r,t)-fir-Ar,t) 



Ar 



(r^Ar) 



Ar2 
24 



(P4,t) + 



r^Ar 



Ar2 
24 



(P5,t) 



Ar^ 
24^ 



60r^^(p3,t) + 60r^^(p3,t) + 15r^^(p3,t) 



Once again, one sees that if /(r, t) is chosen in C(4, 4), the error in this 
difference is bounded. Under this condition, putting all of this information 
together, we expect the error in 

|/(r, t + nAt) - fc{r, t + nAt) \ < KAt, 

where K is bounded and X — > only when both At — > and Ar 0, but 
not when only one of At or Ar go to zero. 

To confirm this, convergence tables have been run on this program. 
/(0, 100) and /(lO, 100) are used as indicators. 

First, the convergence as At — > is investigated. One run is made 
with Ar = 0.100 and At — 0.00125, since we cannot do a run with Ar 
and At infinitely small, and the stability analysis tells us that we must keep 
At < CAr^/^. Call this foo- Then At is allowed to take on a variety of larger 
values and each time we can find the errors £'0 = |/(0, 100) — /oo(0, 100) |, 
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ElO = 1/(10, 100) - /oo(10, 100)1, h = At- 0.00125. Then with any two of 
these one can calculate 

\n{EOa/EOb)/Hha/h) and \n{E10a/E10b)/ Hhjh) (A.23) 

In this problem, this quotient of natural logarithms should be close to 1 since 
in theory the error is KAt. 

The complete set of initial conditions are /(0,t) = 1.0, -Rmax = 100, 



/(0,t) = -0.01. The data is in Table ^ 



In Table |A.1| , skipping the first row and proceeding downward, calculate 



Table [A.2| , where the previous line is the line associated with subscript "a" 



in equation | |A.23| |, and the current line is associated with subscript "b" in 



equation [|A.23|| . It is clear that the quotient of the natural logarithms as in 



A.23II is close to one and gets closer as the size of At decreases, as the theory 
predicts. 

Next, the convergence as Ar -h> is investigated. One run is made 
with Ar = 0.00625 and At = 0.00125, since we cannot do a run with Ar 
and At infinitely small, and consider this to be the value of foo- Then Ar is 
allowed to take on a variety of larger values. Likewise, we calculate the values 



in ||A.23|| , but this time, our theory has no prediction for the value, as the error 



is predicted to decrease as KAt. 

Again, the complete set of initial conditions are /(0,t) = 1.0, -Rmax 
100, /(0,t) = -0.01. The data is in Table [O . 



In Table |A.3|, skipping the first row and proceeding downward, calculate 



Tabel ||A.4]|, where the previous line is the line associated with subscript "a" 



106 



m 



A.23||, and the curent line is associated with subscript "b" in ||A.23||. The 



quotient of natural logarithms is closest to being an integer when the error 
caused by Ar is much greater than the error caused by At. This is as this 
sort of analysis would predict when the error is a sum of a piece that goes to 
zero as Ar and another piece that goes to zero as At 0. 
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Table A.l: 4+1 dimensional model: Convergence data 1. 



Ar 
0.100 
0.100 
0.100 
0.100 
0.100 
0.100 



At 
0.00125 
0.05000 
0.04000 
0.02000 
0.01000 
0.00500 



/(0, 100) 
0.247164210997 
0.247225131669 
0.247212637097 
0.247187644555 
0.247175146883 
0.247168897823 



/(10,100) 
0.245952971169 
0.246014184991 

0.246001630149 
0.245976517632 
0.245963959986 
0.245957680759 



Table A. 2: 4+1 dimensional model: Convergence data 2. 



h 

0.04875 
0.03875 
0.01875 
0.00875 
0.00375 



EO 

0.000060920672 
0.000048426100 

0.000023433558 
0.000010935886 
0.000004686826 



^10 
0.000061213822 
0.000048658980 
0.000023546463 
0.000010988817 
0.000004709590 



In quot. £"0 In quot. ElO 



0.999822247666 

0.999907462503 
0.999972932084 
0.999995539233 



0.999835333093 

0.999894968679 
0.999944166093 
0.999975691565 



Table A. 3: 4+1 dimensional model: Convergence data 3. 



Ar 
0.00250 
0.10000 
0.05000 
0.02500 
0.01250 
0.00625 



At 
0.00125 
0.00125 
0.00125 
0.00125 
0.00125 
0.00125 



/(0,100) 
0.250239554784 
0.247164210997 
0.249463484942 
0.250046460114 
0.250192709322 
0.250229305018 



/(10,100) 
0.248987627424 
0.245952971169 
0.248227707633 
0.248801253159 
0.248943556603 
0.248978377699 
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Table A. 4: 4+1 dimensional model: Convergence data 4. 



h 

0.09750 
0.04750 
0.02250 
0.01000 
0.00375 



EO 

0.003075343787 
0.000776069842 
0.000193094670 
0.000046845462 
0.000010249766 



^10 
0.003034656255 
0.000759919791 
0.000186374265 
0.000044070821 
0.000009249725 



In quot. £"0 In quot. ElO 



1.914735161591 

1.861663705261 
1.746545407713 
1.549300513769 



1.925458076682 

1.880927424171 
1.778154140081 
1.591718519800 



Appendix B 



Stability and convergence of the CP model, 

charge 1 sector 



In this chapter we analyze the stabihty of the equation 



and the associated differencing scheme used in finding numerical solutions. 
This analysis will proceed analogously to that of the 4+1 dimensional model 
outlined in appendix A. Once again / = /(r, t), and r is a radial variable, 
hence r > 0. It will be shown that the stability behavior of the differential 
equation is qualitatively the same as that of the difference equation, and that 
the difference equation converges to the differential equation as Ar —>■ and 
At 0. 



B.l Continuum stability: simplified model 

As in section A.l, the first thing we will address is the stability of the linear 
part of the partial differential equation ||B.1||. This is: 



/ = /" + — (B.2) 

r 



109 



110 



Setting 

f{r,t) = e'^-'e''^', 

and plugging this into |pj.2|] yields the equation: 

2 2 
u) = n . 



Solve this for uj and find that 



uo = ± 



Hence uj always has a negative imaginary part, therefore this equation has a 
growing mode and is not stable, as we saw before with equation ||A.1|| in section 
A.l. 

Once again, we can abstractly consider the equation 

/ = /" + — = 
r 

where Lf is a linear operator with 



and hence that 

£ = r-3/2(r-3/25,r3/2)(r3/2a,r-3/2)r3/2 = -B-'A^AB 

where B = r^^/^ and A = r^^/^S^r^/^. And, just as before A'^A is hermitian 
so it has real spectrum, and since it is essentially a square, it has positive real 
spectrum. Consequently spec(£) = spec(— A'I'A) is real and negative. The 
solutions of the equation 
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consists of sines and cosines in the time variable multiplying the eigenfunctions 
of £, and hence has no growing mode and is strictly stable. 

This is resolved in exactly the same manner as in section A.l. In this 
case we find that under the normal operating conditions of k ^ 1/r we have a 
factor of r~'^l'^ that is not accounted for in the Von Neumann stability analysis 
that had no choice but to appear in the e*"^* portion of the equation. 

When one forcibly puts this factor of r~^/^ into the Von Neumann 
stability analysis via 



and plugs into ||B.2|| , one calculates: 
Hence 

U) = ± , 

2r 

and uj has no negative imaginary part. 



B.2 Continuum stability 

The first thing to address is the stabihty of partial differential equation p.l 



We set 

f{r,t)=fo{r,t)+er-'/'e''^^e'-\ 



and plug into [p.l|| then linearize in e to obtain: 
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3 , 4mr-6 8/^/or Aifofocu 



+ 



4r2 r^ + fi {r^ + f^r + fi 



+ 



r(r2 + /2) r2 + /2 



./I) 



Solving this for cu yields: 

2ifofo 



{r' + Po) 



±K 



1 + 



+ 



3/0' 



14/o/^r 



,,2^2 



K2(r2 + /2)2 /^2(^2 + j2)2 4^2^^2 ^ f2 



2^2 



1/2 



:b.3) 



The imaginary parts are bounded. If we restrict our concern to the 
realm where k is large and -, we actually have that all of the rest of the 



terms in the square root are bounded and small. 



B.3 Discretization scheme 

As in section A. 3, it was found that the main stability problems were generated 
by the linear part of the equation, namely: 

r 

as r ^ 0. As before we use the natural differential operator 
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and we discretize it in the natural way. Letting gAr = r and nAt = t, this 
discretization is 

1 3 



(gAr)" 



q+-] Ar 



f{{q + 1) Ar, nAt) - /(gAr, nAt) 
Ar 



f{qAr, nAt) - f{{q - 1) Ar, nAt) 
Ar 



Ar. 



B.4 Outline of stability analysis 

The stabihty analysis for this equation will be approached in exactly the same 
way as in section A. 4 for the previous model. 

B.5 Stability near zero 

Once again we follow our previous methodologies. A review of section A. 5 
might be in order. One has 

5f^£.^(5f)+L2(5'f), 

with £ji and £2 linear operators. As before £1 and £2 are linear operators 
with £j2 close to a multiple of the identity matrix. The analysis of the eigen- 
values proceeds as before, and once again all hinges on showing that if a is an 
eigenvalue of £1 then a < 0. 

Allow Sf{q,n) to represent 6f{qAr,nAt) and likewise with fo{q,n) = 
/o(gAr, nAt), and compute the following linearized equation, discretized in 
space only: 
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q + - ] Ar^ 



5f{q + l,t)-5f{q,t) ^^ 



Ar 

6fiq,t)-6fiq~l,t) 
Ar 



Ar 



(B.4) 



4/o(g, t)fo{q, t)6f{q, t) ^ 2/o(g, t) {6f{q + 1, t) - Sf{q - 1, t)) 



g2Ar2 + /2(g,t) 



2/o'(g, t) 4/2(g, t)5/(g, t) 



One sees immediately from this that 

4/o/o 



r2 + /2 



cJ 



is a good approximation. 

£-1 is determined by ||B.4 l and the quadratic fit boundary condition at 
the origin, i.e. 

5fiO,t) = pf{Ar,t)-^5fi2Ar,t). 

Let £-1 = [ctij], and /o = fo{r,t) = fo{kAr,t) and obtain the following tridi- 
agonal matrix: 



4/(1/2)= 



2Ar 



3 V Ar2 ' Ar2 + /2 



(3/2)3 (1/2)= 



Ar2 Ar2 
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ai,2 



0'k,k 



0'k,k+l 



_ 1/(1/2)3 
3 V Ar2 



2Ar 



+ 



Ar2 + /2 



+ 



2Ar 



(3/2)^ 

Ar2 Ar2 + 



PAr2 + /2 
+ 1/2)3 (A; -1/2)3 



k^Ar 



(A; + 1/2) 



k^Ar^ 
2k^Ar 



+ 



k^Ar^ A;2 Ar2 + /2. 

One can now use Maple or another program to compute the eigenvalues and 
eigenvectors of this matrix while changing the size of the matrix n and the 
values of /o /o and Ar. 

Since /o is always much less than /o, the contributions from the terms 
with /o are negligible. One expects that rescaling Ar by a factor p should 
result in the eigenvalues changing by approximately 1/p^, since the largest 
terms are multiplied by Ar^. 

Let's see this explicitly. When n = 5,/o = l,/o = — 0.01, Ar = 0.01 
the matrix is: 





-33333.3 33333.3 
4218.8 







■23750.0 19531.2 
5787.2 -21666.7 
6699.5 







15879.4 
-20937.5 
7290.5 







14238.0 
-20600.0 



This matrix has eigenvalues with multiplicity m and eigenvectors given by 



Eigenvalue 


m 




Eij 


j;envector 




-43905.3 


1 


-0.947, 


0.300,- 


-0.105, 


0.038, 


-0.012 


-24319.3 


1 


-1.784, 


-0.482, 


0.399, 


0.109, 


-0.214 


-35118.8 


1 


3.303, 


-0.177,- 


-0.610, 


0.582, 


-0.292 


-12891.4 


1 


2.504, 


1.535, 


0.313,- 


0.387, 


-0.366 


-4052.7 


1 


0.728, 


0.640, 


0.488, 


0.308, 


0.136 
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Now rescale /o = —0.1, which is multiplying by a factor of 10. We 
expect this to have httle or no effect. 



-33333.3 33333.3 

4218.8 -23750.0 19531.2 

5787.2 -21666.7 15879.4 

6699.5 -20937.5 14238.0 

7290.5 -20600.0 



This matrix has eigenvalues with multiplicity m and eigenvectors given by 



Eigenvalue 


m 




Eif 


^envector 




-43905.4 


1 


[-0.947, 


0.300, - 


-0.105, 


0.038, 


-0.012 


-35118.8 


1 


[-3.303, 


0.177, 


0.610, 


-0.582, 


0.292 


-24319.3 


1 


[-1.784,- 


-0.482, 


0.399, 


0.109, 


-0.214 


-12891.4 


1 


[ 2.504, 


1.535, 


0.313, 


-0.387, 


-0.366 


-4052.7 


1 


[ 0.728, 


0.640, 


0.488, 


0.308, 


0.136 



Now rescale Ar = 0.1 so that Ar is 10 times the original with /o = 
1, /o = —0.01. This we expect to scale the eigenvalues by a factor of 100. The 
matrix is: 



-333.1 333.1 

43.0 -237.5 194.5 

59.5 -216.7 

69.8 







157.1 

-209.4 139.6 

76.9 -206.0 



The eigenvalues with multiplicity m and eigenvectors are: 



Eigenvalue 


m 




Eigenvector 




-440.7 


1 


[ 0.945, 


-0.305, 0.110, -.041, 


.0136 


-353.0 


1 


[-3.145, 


0.188, 0.583,-0.576, 


0.302 


-243.4 


1 


[ 1.745, 


0.470,-0.399,-0.110, 


0.226 


-127.6 


1 


[ 2.333, 


1.439, 0.298,-0.376, 


-0.369 


-38.0 


1 


[-0.714, 


-0.633, -0.491, -0.319, 


-0.146 



These are indeed scaled by a factor of 100 from the original. 
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Now rescale so /o = 10 is 10 times as large with Ar = 0.01, /o 
The matrix is: 



-0.01. 



-33333.3 
4218.8 






33333.3 

-23750.0 
5787.0 







19531.2 
-21666.7 
6699.2 








15879.6 
-20937.5 
7290.0 








14238.3 
-20600.0 



The eigenvalues with multiplicity m and eigenvectors are: 



Eigenvalue 
-43905.2 
-35118.5 
-24319.3 
-12891.5 
-4053.0 



m 

1 

1 

1 
1 
1 



Eigenvector 
-0.9471, 0.3004,-0.1054, 0.0381,-0.0119] 
3.3030, -0.1769, -0.6105, 0.5816, -0.2920 ] 

-1.7836,-0.4823, 0.3993, 0.1091,-0.2138] 
2.5037, 1.5354, 0.3128,-0.3867,-0.3657] 
0.7281, 0.6396, 0.4877, 0.3079, 0.1357] 



These are virtually the same as in the original with /o = 1.0. 

The final rescaling is to reset the matrix size n = 10 with Ar 
/o — 1, fo — —0.01. The matrix is: 



-33333. 33333. 
4218.8 -23750. 19531. 



5787.2 -21667. 15879 























6699.5 -20938. 14238. 

7290.5 -20600. 13310. 



7703.3 -20417. 12713. 



8007.5 -20306. 12299. 



8241.0 -20234. 11993. 
8425.8 -20185. 11759. 

8576. -20150. 



.01, 



118 



The eigenvalues with multiplicity m and eigenvectors are: 

Eigenvalue m Eigenvector 
-43977.1 1 

-39218.2 1 

-35858.2 1 



[ 0.9460, 


-0.3021, 


0.1085, 


-0.0423, 


0.0175, 


-0.0075, 


0.0033, 


-0.0015, 


0.0006, 


-0.0002] 






[-2.5044, 


0.4421, 


0.1908, 


-0.3720, 


0.3879, 


-0.3388, 


0.2660, 


-0.1885, 


0.1156, 


-0.0520] 






[-3.5682, 


0.2703, 


0.6032, 


-0.6376, 


0.3843, 


-0.0914, 


-0.1219, 


0.2136, 


-0.1945, 


0.1062] 






[ 6.2021, 


0.4246, 


-1.4984, 


0.7308, 


0.1859, 


-0.5463, 


0.3443, 


0.0549, 


-0.2861, 


0.2250] 






[-2.3840, 


-0.5703, 


0.5619, 


0.0772, 


-0.2884, 


0.0608, 


0.1511, 


-0.1017, 


-0.0604, 


0.0994] 






[-3.8375, 


-1.6147, 


0.4616, 


0.6571, 


-0.1420, 


-0.3737, 


0.0534, 


0.2477, 


-0.01757, 


-0.1788] 






[-4.6228, 


-2.7619, 


-0.4625, 


0.7663, 


0.6223, 


-0.0839, 


-0.4233, 


-0.1824, 


0.1872, 


0.2384] 






[-2.2447, 


-1.6940, 


-0.8658, 


-0.1181, 


0.3015, 


0.3462, 


0.1505, 


-0.0769, 


-0.1808, 


-0.1295] 






[-1.2609, 


-1.1095, 


-0.8494, 


-0.5405, 


-0.2432, 


-0.0072, 


0.1381, 


0.1877, 


0.1592, 


0.0845]} 






[ 0.5060, 


0.4876, 


0.4534, 


0.4063, 


0.3497, 


0.2868, 


0.2214, 


0.1570, 


0.0969, 


0.0439]. 







-31051.4 1 



-25358.7 1 



-19307.4 1 



-13418.0 1 



-8177.3 1 



-4002.9 1 



-1210.9 1 



Once again, all eigenvalues are negative. The modes localized near zero don't 
change. The eigenvalues increase towards zero as n increases because there 
are growing modes of this equation away from r — 0. 

One can conclude that the eigenvalues of this matrix under reasonable 
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initial conditions will always be negative, as required. This differencing scheme 
does not have any instabilities generated at the origin. 



B.6 Stability away from zero: simplified model 

The analysis in this section shall follow that of section A. 6. 

In this section the stability of the difference equation derived from 
||B.1|| with the differencing scheme outlined in section 3.1 under the special 
condition that f{r,t) = fo{r,t) + ge«'"?^»'e«'^"'^* where fo{r,t) = c, a constant, 
will be addressed. 

Explicit inclusion to the factor r~^/^ is omitted, because it does not 
impact this analysis to omit it. 

Plug the previous expression for f{r,t) into [BTT| and linearize in e to 
obtain the following: 



g2Ar2 + /2 ■ 



:b.5) 



J 



As before let x = e"^'^* and 



Reduce equation |p.5|| to 

- (2 + JAt^)x + 1 = 0. 
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This is exactly as in section A. 6. As before there are two logarithmic solutions 
for CO given by 



UJ — 



-i In 1 + — — ± VJAt\ 1 H — 



At 



As in section A. 6, the biggest concern is determining exactly how big 
JAt^ is. Let kAt = 6, then reduce exponentials to sines and cosines appropri- 
ately, and expand the factors of g ± 1/2, and simplify using such information 
as q — r/Ar, to obtain 



2cos^-2\ 3m /sin ^\ 3fi;2ArV2 cos^ - 2 



02 J r- \ 9 J Ar^ \ 9^ 
ikAt"^ ( sin^\ ^ir^K ( sin^' 



4r3 V ^ y r2 + /2 V ^ 
As At — > one of these terms go to zero but the following do not: 

2 I 2cos^ - 2\ 3m /sin^\ S/t^Ar^ /2cos^ - 2\ ^ir'^K f sm9 



Q2 J ' ^ y J ' 4^2 \^ 02 J j.2j^j2\^Q 

As before, the trigonometric parts of these terms is strictly bounded, and we 
are working in the case where r is bounded away from zero. Our only concern 
is with the size of ac, but we know that nAr < tt and At <C Ar, hence we 
can make At < cAr and with c chosen sufficiently small, we can make nAt 
as small as we like. In particular we may use the Taylor approximation to the 
square root and then apply the Taylor approximation to the logarithm. Just 
as in section A. 6, with some cancellations, one obtains: 



—i In 

UJ fti 



± ^/jAt i 1 + + 0{J^At^] 
At 
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—I 



±\fjAt + O ({\fJAt) 



At 



+ (yi^At^^ 



Recall, as Ar —>■ all of the terms in J go to except for 



2 cos 6' - 2 



So one has 



3iK / sin^ 



S/t^Ar^ /2cos0 - 2 



4 J" 2 



^2 



Air'^K, ( sin^ 



r2 + /2 



2cos^ - 2 

^2 



+ 



3m /sin^\ S/c^Ar^ /2cos^ — 2 



+ 



4^2 



^2 



r2 + /2 



+ 0(Ar) 



Since is large while Ar is small, approximating this square root yields: 

^2 



2 - 2cos^ 
iAr^ / sin^ 



1 + 



3i 
2Kr 

^2 



sm I 



2 - 2cos( 



3AH 



2ir' 



/ sin ^ 



2-2cos^y /«(r2 + /2) V ^ yV2-2cos^ 



^2 



+0 



As before, the leading order term of \fj is 0{k). We can make the correction 
0(\/j At^) as small as we hke by requiring that At < cAr^/^ ^^d choosing 
the factor c appropriately. Given these choices, one has: 



LiJ 



± 



2 - 2cos^ 



3i / sin^ 
+ — 



2r\Q 



iAr^ ( sin^ 



^2 



8r2 



2 - 2cos6' 
^2 



+ 



3KAr2 



8^2 



2ir2 



2-2cos^y (r^ + /o) 



sin^ 



^2 



2 - 2cos^ 
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It is also clear as in section A. 6 that if one permitted r ^ the dominant 
piece would be: 



/2-2COS0 3i fsm0\ ( 0^ 
UJ ^ ±fi;\/ =F — 



02 ' 2r V / V2 - 2cos0 
and as ^ this goes to 

2r 

or exactly what we obtained in section B.l that was so often "corrected" by 
the addition of the factor r^^/2_ 



B.7 The complications 



We must now deal with the reality that the assumption /o = c is not in general 
true. We shall, as in section A. 7, analyze the difference between the solution 
for u when /o = c, called c^o, and the solution for u when /o ^ c, called ui. 

To find ujo we solved ||B.5|| an equation which can be characterized as 

g{iu) = old(K). 



It was: 



oAt 



At2 



(.+ 0'(e---l)-(g-i)'(l-, 
g^Ar2 



-iKAr\ 



g2Ar2 + /2 



If we do not assume that fo = c then the same linearization in e of [[B.l|l with 



/(gAr, nAt) = /o(gAr, nAt) + ^g^^Argit^nAt yjgjj^g ^j^g following equation: 

3 / 1 \ 3 



oAt 



±e_ 

At^ 



oAt 



(1- 



-iKAr\ 
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2At(r2 + /2) 2Ar(r2 + /2) 

2 (/o^ - a _ mitiii' 

+ Po {r^ + /o )^ 



Characterize this equation as 

g{w) — old(«;) + newi(«;) + new2(a;) 



with 



"^^^^^^ - 2Ar(r2 + /2) (r^ + /2)2 



^2 + f2 (^2 + J2)2 



and 



(B.7) 



I X , JoJole - e J 

"^"^("^ = + 2At(r2 + /2) ^^-S) 

Here /o is substituted for the finite difference 

/o(gAr, (n + 1) At) - /o(gAr, (n - 1) At) 
2At 

and hkewise with /q. 

Prom here we proceed as in section A. 7. We have 

, , 2 costuAt — 2 ,^ 
9i^) = ^2 ' (B-9) 

•^0 = 5'~^(old(K)), 

and 

a;i = g'~"^(old(«;) + newi(«;) + new2(a;i)). 



124 



Let 



then 



ujs = g ^(old(K) + s(newi(K) + new2(co's))) 



\UJi - UJo\ 



ds 



ds 



< \[{9 )'](newi(K) + new2(u;J)| (is 
Jo 

< max (K^f^"*^)')!) max (|newi(K) + new2(co's)|) (B.IO) 



So finding a bound on uj in terms of n will allow one to find a bound for 
\uJi — ujq\ in terms of k. 

We start by showing that 

T(uj) = g^^{old{K) + newi(K) + new2(co')) 

is a contraction mapping on the ball Icj— k| < \k,\/2 = B{k,, \k,\/2). Then we will 
know that there exists a solution for u G B{k,, |k|/2) and hence \uj\ < 3|/t|/2. 



As before using |B.9 



-2sin(cjAt) 
At ■ 



We've already stated one needs to choose At < cAr^/^ in section B.6, so 
specifically, one can make At < Ar/150. Then since |KAr| < tt and |a;| < 
3|k|/2 this implies that IcijAtl < vr/lOO and one has 

sin(ti;At) 



\g'{uj)\ = \2uj\ 



> 1.98L' . 
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If LUi and LJ2 are in B{k,, k,/2), assign yi and y2 as follows and find that 



yi — old(fi;) + newi(K) + new2(a;i) 



Letting, 



calculate 



y2 = old(K) + newi(fi;) + new2(cc;2) 
y2-y\ = new2(a;2) - new2(a;i). 



y{s) =1/1 + 5(2/2-1/1), 



\T{U2)-TM\ = 



\9-\y2) - 9-\y,)\ 



{9 ) {y{s)){y2 - yi)ds 



< [ max\{g ^)'{y{s))\\{y2-yi)\ds 
Jo 

< max |(r7(2/(s))| 1(2/2 -t/i)|. (B.ll) 



Simple calculus yields: 

(r^)'(^(^)) 



9'H' 



By choice of yi and 2/2 and 001,002 G B{k, k/2), 



max I (5 ^)'(2/(s))| = — 



< 



.98\u\ - 1.98\k\ 

Analyzing since |a;i|, \u!2\ < 3|k|/2 and |a;iAi|, |a;2Ai| < vr/lOO, 

/o/o sin(a;2 At) isin(a;iAt) 



(B.12) 



|new2(a;2) - new2(a;i)| 



r2 + /2 



At 



At 
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< 




r2 + /2 

Now plugging |B.12| and [[B.13[| into |B.11|], obtain: 



sin{uj2At) sin(u;iAt)~ 



At At 
cos{sAt)ds 



\UJ2 — UJl\ 



(B.13) 



\T{uj2)-T{uj,)\ = \g-\y2)-g-\yi)\< 



1.98k| 



/o/o 



r2 + /2 



The quantity 



/o/o 



/o + r2 

defines the time scale for the problem. In order for this entire analysis to make 
sense one expects 

K ^ time scale and k ^ length scale. 

Therefore, this is sufficient to show that T{w) is a contraction map from 
B{k,k/2) to itself. So by the Contraction Mapping Principle, one concludes 
that there exists a fixed point of T in B{h,k,/2). Clearly this argument 
holds just as well for — T and B{—k,,k,/2), so there is a fixed point for —T 
in B{—K, k/2). Since the left hand side of |[A.7|| is quadratic in e*"^^*, there are 
two solutions for u and these are they. Lastly, we finish estimating in equation 
|B.10| |. We already have a perfectly good estimate for {g'^)' from |B.12| . We 
merely need to estimate 



newi(fi;) + new2{uj) \ . 
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Replacing the exponentials with their trigonometric forms in [p.7|| and ||B.8 
to get: 

ifofo f sin(KAr 



newil Kl 



(r2 + 



cl2\ 



r2 + /2 



and 



new2(u; ) 



(r2 + 
/o/o /isin(c^At) 



r2 + /2 



At 



Recall also that K^l/rso/tr^l. Now if 6' is real, 

Isin^l < e, 

and if 9 is complex with \9\ sufficiently small, as it would be if ^ = cuAt, then 

Isin^l < l.Oll^l, 



one has 

|newi(K) + new2(ci;) 



< 



/n/o 



r2 + /2 



k + 



+ 



,2 + J2)2 

fofo 



+ 



2(/o^ - 



'2^ 



+ 



r2 + /2 



+ Po 



|1.01u;| 



(r2 + fif 

By my previous arguments about n and the length and time scales, 
the three central terms are much much smaller than the others. Plugging 
\^\ < 3|k|/2 and this and ||B.12|| into ||B.10|| one concludes that under these 
assumptions 



\Wi-Wq\ < 



1.98|/t| 



/o/o 



r2 + /2 



\K\ 



(r2 + 



+ 



2(/o^ - /^^) 



+ 



+ 



4/o^(/o' - /^^) 



(r2 + 



/o/o 



r2 + /2 



3.03|k| 
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Hence, \uji — ujqI is bounded, which imphes in turn that the imaginary part 
of uji is bounded since the imaginary part of ujq was. Once again the growing 
modes are bounded. 



Further, if we compare with equation] B.3 |, we see analogues between 
the terms in these two equations. 



1.98 



/o/o 



and 



/2 + r2 

3.03 

L98 



+ 



r' + fi {r^ + Po? 



/o/o 



/^ 



-2^/o/o 



B.8 Convergence estimates 

This section contains an analysis of the convergence of the differencing scheme 
for the equation 

2/ 



/2 



(B.14) 



We will show that as Ar —>■ and At that the solution found for f{r,t) 
converges to an actual solution of the partial differential equation. As in 
this entire chapter, the arguments will essentially mimic those given before in 
section A. 7. 



Substituting the appropriate differences into |B.14| | and forward inte- 
grating the equation one can solve for fc{r,t + At), where fc is used instead 
of f{r,t) to indicate it is a calculated value (at (r, t)). One obtains 



/c(r,t + At) 



2/c(r,t)-/e(r,t- At) 
3/^ 4r/^ 



+At^ 



/:+ 



+ 



2/c 
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Here we also assume the derivatives are represented by the appropriate differ- 
encs as outhned in previous sections. Approximate /c(r, t + At) for use in 
by using first fc{r,t + At) = 2/c(r, t) — fc{r,t — At) and then iterating the 
solution found for /c(r, t + At) six times, finding a new one from solving with 
the previous value. Using a greater number of iterations on /^(r, t + At) does 
not change the answer to the numerical precision on the computer, so this is 
not considered a source of error. 

The error comes from the discretizations. If we let At^Ci be the error in 
the difference for / and C2 be the accumulated error in the differences for the 
left hand side of |[B.14|| , then one may complete the exact same computations 



as in section A. 8 to find the error in n steps is bounded by 

|/(r,t + nAt) -/,(r,t + nAt)| < nAtWCi\ + \C2\) 

+ AtC3(|Ci| + |C2|) 

where C3 is a constant. 

Taylor's theorem will be applied extensively to the terms of equation 
B.14II to show that Ci and C2 are constants. We see 



■■ /(r,t + At) + /(r,t-At)-2/(r,t) At'd'f 
nr,t) = — 12- a^^^'^^) 

4r/'(r,t) 2(/(r + Ar,t) -/(r- Ar,t)) 



/2(r,t)+r2 Ar(/2(r,t) +r2) 

4r fAr^d^f 



2f{r,t)p{r,t) ^ f{r,t){f{r,t + At)-f{r,t-At)f 
P{r,t)+r^ 2At2(/2(r,t) + r2) 
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2/(r,t)r(r,t) 
y2 _|_ J.2 



2f{r,t){f{r,t + At)-f{r,t-At)) 



At(/2(r,t) +r2) 



X 



At2 d^f 



r, m) 



6 df-^ 
2/(r,t) /At^^V, 



/2(r,t) +r2 V 6 (9t3 

/(r,t)(/(r + Ar,t)-/(r-Ar,t)f 
2Ar2(/2(r, + r2) 

2/(r,t)(/(r + Ar,t)-/(r-Ar,i)) 



Ar{p{r, t)+r^) 



X 



Ar^ 



6 

2/(r,i) fAr^d^f 



/2(r, + r2 V 6 ^r^ 



The first of these equations, for example, says that if 



(r, Tji) 



is bounded then the error Ci — > as At'^ — > 0. If f{r,t) is chosen in C(4, 4), 
the space of functions with five bounded space/time derivatives, all of these 
differences have bounded errors. 



Now we must attend to the discretization error for 



r(r,t) + 



3/'(r,t) 



r + 



- r - 



Ary f f{r + Ar,t) - f{r,t) 
Ar 

ArV ff {r,t)-f{r-Ar,t) 
Ar 



(B.15) 



(r^Ar). 



If we let h(r) — r f'(r) we have: 



h'{r) = 3r^f{r) + r^f"{r) = /"(r) + 



3/'(0 
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h"{r) = 6r/'(r) + 6rV"(r) +r7"'(r) 



h"'{r) = 6/'(r) + 18r/"(r) + 9r^f"'{r) + r''f^\r). 
Using the appropriate Taylor Series, one can calculate that 

„3 



h'{r) 



h{r + ^) - h{r - ^) Ar- 



Ar 



24 



Now if we allow 



Pi e 



Ar Ar 

r , r H 

2 ' 2 



and substitute, to obtain 

3/'(r) 



fir) + 



Ar 



Ar^ 



(6/'(pi) + 18rr(pi) + 9rV'"(pi) + r^^^Hpi)) (B.16) 



Now calculate central difference approximations for f'{r + Ar/2) and 



/'(r — Ar/2) and substitute into ||B.16|| . Let p2 G [r, r+Ar] and pa G [r — Ar, r]. 
The resultant equation is 
3f{r,t) 



f"{r,t) + 



r 



r + 



ArV ff{r + Ar,t) - f{r,t) 



(r + Al) 



2 y V Ar 

3 



ArV //(r,t) - /(r - Ar,t) 



Ar 



r^Ar 



Ar^ 
24 (9r3 



(P2,t) + 



(r^Ar 



Ar2 
24 9r3 



(P3,i) 



Ar^ 
'2^ 



(6/'(pi) + 18rr(pi) + 9rV'"(pi) + r^/^^HPi)) • 
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Once again one sees that if f{r,t) is chosen in C(4, 4), the error in this differ- 
ence is bounded. Under this condition, putting all of this information together, 
we expect the error in 

\f{r,t + nAt) - fc{r, t + nAt) \ < KAt 

where K is bounded and K Q only when both At and Ar 0, but 
not when only one of At or Ar go to zero. 

To confirm this, convergence tables have been created for this program. 
/(0,60) and /(10,60) are used as indicators. 

First, the convergence as At ^ is investigated. One run is made 
with Ar = 0.100 and At = 0.001, since we cannot do a run with Ar and At 
infinitely small, and the stability analysis tells us that we must keep At < 
CAr3/2. Qg^u ^j^jg r^Yien At is allowed to take on a variety of larger values 
and each time we can find the errors = |/(0,60) — /oo(0,60)|, ElO = 
1/(10, 60) - /oo(10, 60)1, /i = At - 0.001. Then with any two of these one can 
calculate 

\ji{EOjEOb)/Hha/h) and HElOjElOb)/ Hhjh) (B.17) 

In this problem, this quotient of natural logarithms should be close to 1 since 
in theory the error is KAt. 

The complete set of initial conditions are /(0,t) = 1.0, -Rmax = 100, 
/(0,t) = -0.01. The data is in Table PTT. 



In Table p.l| , skipping the first row and proceeding downward, calculate 



Table p.2| , where the previous line is the line associated with subscript "a" in 
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equation p.l7|| , and the curent line is associated with subscript "b" in equation 

is 



B.17|| . It is clear that the quotient of the natural logarithms as in [|B.17 



close to one and gets closer as the size of At decreases, as the theory predicts. 

Next, the convergence as Ar — > is investigated. One run is made with 
Ar = 0.01 and At = 0.001, since we cannot do a run with Ar and At infinitely 
small, and consider this to be the value of foo- Then Ar is allowed to take on 
a variety of larger values. Likewise, we calculate the values in [|B.17|] , but this 
time, our theory has no prediction for the value, as the error is predicted to 
decrease as KAt. 

Again, the complete set of initial conditions are /(O, t) = 1.0, -Rmax = 



100, /(0,t) = -0.01. The data is in Table ^ 



In Table [B.3| , skipping the first row and proceeding downward, calculate 



Table B.4, where the previous line is the line associated with subscript "a" in 



equation |p.l7|| , and the curent line is associated with subscript "b" in equation 
||B.17|| . The quotient of natural logarithms is closest to being an integer when 
the error caused by Ar is much greater than the error caused by At. This is 
as this sort of analysis would predict when the error is a sum of a piece that 
goes to zero as Ar and another piece that goes to zero as At — * 0. 
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Table B.l: CP^ model, charge 1 sector: Convergence data 1. 



Ar 


At 


/(0,60) 




/(10,60) 


0.1 


0.001 


0.433979230592 


0, 


.431548378650 


0.1 


0.0025 


0.434046902945 


0, 


.431616258041 


0.1 


0.005 


0.434159692320 


0, 


.431729392396 


0.1 


0.0075 


0.434272484358 


0, 


,431842529337 


0.1 


0.0125 


0.434498076252 


0. 


.432068810954 



Table B.2: CP^ model, charge 1 sector: Convergence data 2. 

In quot. EO In quot. ElO 



h EO 

0.0015 0.0000676724 

0.0040 0.0001804617 

0.0065 0.0002932538 

0.0115 0.0005188457 



ElO 
0.0000678794 

0.0001810137 
0.0002941507 
0.0005204323 



1.000011111 
1.000024846 
1.000040720 



1.000011077 
1.000024070 
0.994688335 



Table B.3: CP^ model, charge 1 sector: Convergence data 3. 



Ar At /(0,60) 

0.010 0.001 0.434213647771 

0.020 0.001 0.434206452051 

0.025 0.001 0.434201061431 

0.040 0.001 0.434177743230 

0.050 0.001 0.434156269734 

0.100 0.001 0.433979230592 

0.200 0.001 0.433312345450 



/(10,60) 
.431709360602 
.431707448176 
.431704981289 
431691546698 
431677574831 
.431548378650 
431017145989 
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Table B.4: CP^ model, charge 1 sector: Convergence data 4. 

In quot. EO In quot. ElO 



h EO 

0.010 0.0000071957 

0.015 0.0000125863 

0.030 0.0000359045 

0.040 0.0000573780 

0.090 0.0002344172 

0.190 0.0009013023 



ElO 
0.0000019124 

0.0000043793 
0.0000178139 
0.0000317858 
0.0001609820 
0.0006922146 



1.37897263 
1.512310436 
1.629570850 
1.735588896 
1.802345190 



2.043406152 
2.024231211 
2.012779570 
2.000508645 
1.952054808 



Bibliography 



[1] M. F. Atiyah. Geometry of Yang-Mills Fields. Accademia Nazionale Dei 
Lincei Scuola Normale Superiore, 1979. 

[2] Robert Leese. Low-energy scattering of solitons in the CP^ model. Nu- 
clear Physics B, 344:33-72, 1990. 

[3] B. Piette and W. J. Zakrzewski. Shrinking of sohtons in the (2+1)- 
dimensional S"^ sigma model. Nonlinearity, 9:897-910, 1996. 

[4] Wojciech J. Zakrzewski Robert A. Leese, Michel Peyrard. Soliton stability 
in the 0(3)(T model in (2-1-1) dimensions. Nonlinearity, 3:387-412, 1990. 

[5] J. M. Speight. Low-energy dynamics of a CP^ lump on the sphere. J. 
Math. Phys., 36:796, 1995. 

[6] D. Stuart. The geodesic approximation for the Yang-Mills-Higgs equa- 
tions. Commun. Math. Phys., 166:149-190, 1994. 

[7] R. S. Ward. Slowly-moving lumps in the CP^ model in (2-1-1) dimensions. 

Physics Letters, 158B(5):424-428, 1985. 

[8] Wojciech J. Zakrzewski. Soliton-like scattering in the 0{3)a model in 
(2-1-1) dimensions. Nonlinearity, 4:429-475, 1991. 



136 



Vita 



Jean-Marie Linhart was born in Melrose Park, Illinois on June 2, 1969, the 
daughter of Harriet Roeters Linhart and Carl George Linhart. After complet- 
ing her work at Oak Park and River Forest High School, Oak Park, lUinois, 
in 1987, she entered The University of Chicago in Chicago, Illinois. She re- 
ceived the degree of Bachelor of Science in Mathematics from the University 
of Chicago in June 1990. She entered the Graduate School of The University 
of Texas at Austin in August of 1990. During the summer of 1991 she was 
employed at Lawrence Livermore National Laboratories in the Global Climate 
Research group. While working on her Master's degree, she worked as a Teach- 
ing Assistant in the Department of Mathematics at the University of Texas 
at Austin. She received the degree of Master of Arts in Mathematics from 
the University of Texas in May 1993. She then worked at Lawrence Livermore 
National Laboratory in the Global Chmate Research group until June of 1994. 
Prom July of 1994 to October of 1995, she was employed as a Mathematics 
Instructor at ITT Technical Institute in Austin Texas. In August of 1995, she 
returned to the Graduate School of The University of Texas. Prom August 
1995 to August 1997, she was employed in the Department of Mathematics 
at The University of Texas at Austin as a Teaching Assistant and Assistant 
Instructor. During the summer of 1996 she was employed by Schlumberger 
Well Services in the Formation Evaluation department. She returned to em- 



137 



138 

ployment by Schlumberger Well Services in August of 1997. 



Permanent address: 2207 S. Fifth Street #111 
Austin, TX 78704 



This dissertation was typeset with I^T[t;X|3 by the author. 



^I^T^T^X is a document preparation system developed by Leslie Lamport as a special 
version of Donald Knuth's Tj^X Program. 



